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PREFACE 

The  following  report  has  been  prepared  in  accordance  with  the  terms 
of  Contract  No.  DA -36 -03^ -OED- 1614-6  and  constitutes  the  Final  Eeportccall^d 
for  under  the  terms  of  that  contract. 

Said  Contract  between  the  Institute  for  Advanced  Study  and  the  Depart- 
ment of  the  Army  was  entered  on  June  2,  195^ »  ^or   "further  development  of 
principals  and  methods  for  operation  and  maintenance  of  very  high  speed 
digital  electric  computer  devices",  i.e.  for  continuation  of  our  work  under 
Contract  No.  DA-36-031+-OED-1330,  which  terminated  on  June  30,  1951+,  and  for 
which  a  final  report  was  submitted,  dated  December  195^. 

Contract  No.  DA-36-03i|-CiRD-l61+6,  together  with  its  supplemental  agree- 
ments No.  1  and  No.  2,  was  in  effect  from  July  1,  195^  through  October  31, 
1956  and  was  the  only  contract  supporting  maintenance  and  operation  of  our 
Computer  during  that  period,  except  for  the  last  7  months  during  which  a 
small  part  of  the  operation  hag  been  paid  for  under  Contract  Nonr  1358- (O**-). 
This  new  contract  between  the  Institute  for  Advanced  Study  and  the  Office 
of  Naval  Research  has  taken  over  the  full  burden  of  machine  maintenance  and 
operation  on  January  1,  1957« 

This  Final  Eeport  is  divided  into  two  parts: 

PART  I  covers  the  engineering  work  carried  out  from  July  1,  195U- 
through  December  31,  1956  under  the  terms  of  Contract  No  DA-36-03'i--C!RD-l6l4-6. 

PART  II  lists  a  number  of  problems  for  which,  during  the  same  30 
months,  numerical  results  have  been  obtained  with  the  help  of  our  Computer. 
This  work  was  supported  by  four  contracts; 

(1)  Contract  DA-36-03'+-aRD-l6U-6  and  Nonr  1358-(Ol4-)  for  machine 
operation. 

(2)  Contract  N7-onr-388  for  the  development  of  methods  for  high-speed 
automatic  computing  (from  July  1,   I95J+  through  December  31,  195^.) 

(3)  Contract  Nonr  1358- (03)  for  the  development  of  methods  for  high- 
speed automatic  computing  (since  January  1,  1955). 

(U)  Contract  Nonr  1358-(02)  for  all  mathematical  and  coding  work 
connected  with  meteorological  research. 

During  the  whole  period  our  Computer  was  used  for  scientific  computa- 
tions exclusively  and  no  charge  was  made  to  any  other  organization  or  contract 
for  coding  or  machine  operation . 


Hans  J.  Maehly 
Acting  Project  Director 
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I.   GENEBAL  PURPOSE   ROUTINES 


10.  General  Bemarks 

It  was  decided  early  in  1956  that  no  new  machine  should  be  built  at 
the  Institute  for  Advanced  Study;  that  aost  of  the  engineering  staff  would, 
therefore,  leave  to  pursue  development  work  at  other  places;  and  that  the 
Electronic  Computer  should  be  transfonssed  from  an  experimental  project  into 
a  tool  for  the  solution  of  the  many  computational  problems  arising  in  the 
scientific  community  of  Princeton.  This  transformation  will  be  no  minor 
task,  as  no  attempt  has  been  made,  during  the  construction  and  completion 
of  our  equipment,  to  provide  for  easy  communication  between  the  computing 
machine  and  the  coder  or  mathematician;  and  it  is  made  even  more  difficult 
by  the  requirement  that,  lacking  funds  for  a  staff  of  professional  coders, 
coding  procedures  should  be  simplified  so  that  most  problems  can  be  coded 
by  their  originators. 

After  1  July  1956,  the  development  of  a  consistent  system  of  general 
purpose  routines  has,  therefore,  been  given  priority  over  all  new  work,  which 
simply  means  that  no  major  help  could  be  rendered  to  anybody  presenting 
problems  after  this  date.  We  felt  obliged,  however,  to  finish  those  problems 
whose  coding  was  already  under  way.  This,  together  with  the  engineering 
changes  and  difficulties  described  in  voltune  I  of  this  report,  has  delayed 
the  development  of  our  general  purpose  routines  much  more  than  originally 
anticipated » 

Fortunately,  the  floating-point  interpretive  routine  (FLINT)  went  al- 
ready into  operation  late  in  1955  and,  despite  its  admitted  drawbacks,  has 
largely  filled  the  gap  between  the  restriction  of  our  coding  services  and 
the  completion  of  other  service  routines.  Besides  its  operational  applica- 
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tions,  the  modified  interpreter  (chapter  11)  has  yielded  valuable  sugges- 
tions for  an  improved  version  which  is  still  being  planned.  At  this  point, 

however,  emphasis  is  being  put  on  the  construction  of  subroutines,  up  to 

*) 
a  symbolic  assembler  ,  to  facilitate  coding  in  direct,  fixed-point  machine 

language  (chapter  13) • 

Our  staff  has  been  too  small  to  engage  in  any  ambitious  projects  involv- 
ing new  concepts.  What  is  now  most  needed,  are  standard  utility  routines, 
which,  unfortunately,  we  cannot  copy  from  any  other  place  since  no  machine 
resembles  ours  closely  enough.  A  description  of  our  new  service  routines 
is  included  in  this  report  as  some  of  the  details  may  be  new,  or  of  interest 
where  the  construction  of  similar  codes  or  machines  is  being  contemplated. 


*)  To  be  described  in  Final  Report  for  Contract  Nonr-1358-(0^) . 
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11.  FLINT  =  FLOATING  TOim   INTERPRETIVE  SUBBOUTINE 
Originator  and  Coder:  Hans  J.  Maehly 

11.0  Purpose  and  Applicationa 

During  the  construction  of  the  IAS  Computer,  emphasis  was  on  speed 
and  flexibility  rather  than  on  ease  of  coding.   It  was  thus  possible  to 
carry  out  problems  which  had  been  too  long  or  too  unusual  in  their  kind 
for  earlier  computers.   It  must  be  admitted,  on  the  other  hand,  that 
routine  problems  for  which  many  other  computers  would  have  been  quite 
adequate  took  an  undue  amount  of  coding  time,  while  the  saving  in  runn- 
ing time  (e.g.  2  rather  than  10  minutes)  was  insignificant.  Such  small 
problems  did  arise  in  the  Princeton  scientific  community  and  it  would 
have  been  impractical  not  to  solve  them  here  as  long  as  computer  time 
was  available . 

An  auxiliary  routine  has  therefore  been  written  which,  as  far  as 
its  user  is  concerned,  transforms  our  machine  into  a  slower,  less  sophis- 
ticated instrument  for  which  coding  is  much  simpler.  Meanwhile,  experience 
has  sho\m  that  the  effort  to  write  this  routine  (approximately  one  half 
man-year)  was  a  well  worth-while  investment.  FLINT  has  been  used  not  only 
by  an  increasing  number  of  occasional  users  to  code  their  own  problems 
but  also  by  the  permanent  staff  of  our  organization  for  several  auxiliary 
investigations  and  even  for  parts  of  large  problems. 

11.1  Basic  Structure 

The  basic  features   of  FLINT  are   indicated  by  its  name:      floating-point 
interpretive  subroutine. 

Each  floating-point  number  is  stored  in  one  word  =  ko  binary  bits; 
the  first  31  bits  are  reserved  for  the  "mantissa"  r  and  the  last  9  bits 
for  the   "exponent"©    .     Their  respective  ranges  are: 

|<lr|  <    1 
-256  <    ^     <  256  ' 


*'       Actually  256+9  is   stored  in  the   last  9  bits   in  order  to  prevent  carry 
into  the  mantissa  part. 
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In  all  operations  (such  as  +,  -,  x,   -^)  this  word  is  treated  as  a  number  x 

X  =  r.2f  , 
hence      2"^^^  <  |xl  <  2"^^^^ 
or  appr.  lO"'^'''  <  |x  |  <  5.I0''' 

This  range  is  wide  enough  to  avoid  scaling  problems  for  all  but  the 
most  extreme  cases  (such  as  combinatorial  problems  with  more  than  58.'  or 
some  astrophysical  computations).  The  remaining  31  binary  digits  for  the 
mantissa  still  yield  a  relative  accuracy  of  9  decimal  digits. 

An  interpretive  routine  is,  by  definition,  a  code  that  "translates" 
orders  given  in  a  new  "language"  into  ordinary  "machine  language".  Thus, 
every  FLINT  order  is  picked  up  by  FLINT  as  soon  as  the  previous  order 
has  been  executed,  is  then  interpreted  and  causes  a  series  of  machine 
operations  which  provide  the  changes  required  by  the  given  FLINT  order. 
For  example:  a  FLINT  multiplication  order  will  cause  FLINT  to  multiply 
the  mantissas  and  to  add  the  exponents  of  the  two  respective  floating- 
point numbers.  Thus  the  machine  plus  FLINT  will  act  like  a  new  machine 
though  no  physical  changes  have  been  made  for  that  purpose.  We  shall, 
therefore,  speak  of  FLINT  as  if  it  were  a  virtual  machine  rather  than 
an  auxiliary  code. 

The  length  of  an  order  in  FLINT  has  been  cut  down  to  10  bits  (as 
compared  with  20  for  "direct  coding"),  vhich  are  divided  into  two  pentades-- 
00000  through  11111,  i.e.  0  through  31.  The  first  pentade  defines  the 
operation,  the  second  one  gives  the  address  of  the  operand.  Each  code 
written  for  the  use  with  FLINT  must  be  segmented  into  blocks  of  less 
than  30  words  (120  orders);  all  blocks  will  be  transferred  to  the  magnetic 
drum  during  the  read -in  and  only  one  at  a  time  will,  for  execution,  be 
called  into  the  Williams  Memory  where  line  30  has  been  set  aside  for  this 
purpose.  These  blocks  are,  therefore,  often  referred  to  as  "lines"  of 
FLINT  code. 

Another  transfer  from  the  drum  is  necessary  to  invoke  either  of  the 
two  subroutines  which  provide  for  the  conversion  of  floating  point  decimal 
numbers,  punched  on  and  read  in  from  IBM  cards,  to  floating  point  binary 
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form;  or  for  the  reconversion  from  binary  to  decimal  for  punching  out. 

All  of  these  operations  --  transfer  to  another  line  of  FLINT  code, 
conversion,  and  reconversion  --  are  initiated  by  a  single  order  and 
include  automatic  checking  and,  for  the  latter  two  orders,  the  necessary 
IBM  operations  (read-in  /  punch-out) . 

As  an  aid  in  finding  errors  in  one's  program,  FLIM"  is  available  in 
a  Tracing  mode  in  which  the  contents  of  the  X  and  A  accumulators  and  the 
Y  register,  defined  in  the  following  paragraphs,  are  punched  out,  together 
with  the  address  of  the  instruction  word  to  which  control  is  about  to  be 
transferred.  This  punching  takes  place  at  the  end  of  every  instruction 
word,  i.e.  after  about  every  four  orders.  A  more  flexible  Tracer  is  being 
prepared,  but  the  current  one  has  been  found  to  be  quite  valuable. 

11.2  "Direct"  Floating  Arithmetic  Orders 

FLINT  is  a  single -address  system.  Most  of  the  arithmetic  that  FLINT 
performs  on  the  customer's  data  will  be  done  in  an  accumulator  which  we 
refer  to  by  the  letter  X.  When  not  being  processed  by  the  accumulator, 
data  are  stored  in  floating-point  registers,  m.  The  contents  of  any 
register  are  denoted  by  (  ) ,  and  in  the  case  of  the  accumulator  X,  they 
may  also  be  referred  to  as  x.  The  fundamental  arithmetic  orders  are 
Read  (into  the  accumulator).  Add,  Subtract,  Multiply,  Divide,  and  Store 
(in  floating-point  storage  registers).  The  orders  which  cause  these 
operations  are  shown  below; 

(10, m)         (m)  ^>  X 

(11, m)      X  +  (m)  >  X 

(12, m)      X  -  (m)  >  X 

(13, ra)      X  .  (m)  >  X 

ilk,m)  X  /  (m)  >  X 

(15,m)  X  >  m,  non-clear  X 

In  each  of  these  orders  we  may  refer  by  the  number  m  to  any  of  some  thirty 
data  registers.   (The  binary  character  of  the  IAS  machine  shows  through  in 
the  FLINT  code  by  the  great  importance  attached  to  the  number  32.  The  order 
codes  may  take  on  the  values  from  0  to  31*  Groups  of  addresses  are  treated 
systematically  except  possibly  for  the  numbers  at  each  end  of  the  sequence. 
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i.e.  the  addresses  0  to  31,  and  occasionally  30.) 

As  an  example  of  these  orders,  the  order  (11,17)  would  cause  the  float- 
ing-point number  in  memory  register  17  to  be  added  to  whatever  was  current- 
ly standing  in  X  and,  of  course,  it  would  not  disturb  the  contents  of  regist- 
er 17.  Similarly,  the  order  (15,1^)  would  cause  the  contents  of  the  accu- 
mulator X  to  be  stored  in  register  ik   and  would  not  disturb  the  contents 
of  the  accumulator. 

Although  the  arithmetic  orders  already  given  will  suffice  to  carry 
out  any  operations  that  may  be  needed,  it  is  convenient  to  have  a  few 
more.  The  orders  (l6,m)  through  (19, m)  cause  FLINT  to  execute  the  instruc- 
tions shown  below: 

(l6,m)       (m)   >  Y 

(17,m)   X  +  (m).y  >  X 

(18, m)   X  -  (m)  .y  >  X 

(19,m)    (m)  /  X   ^>  X 

The  order  (l6,m)  is  a"Fetch" order  that  brings  a  floating-point  number 
from  register  m  into  a  special  register  designated  by  the  letter  Y  whose 
contents  we  denote  --  in  analogy  to  the  X  accumulator  --  by  y.  This  register 
is,  in  fact,  register  0  --  so  the  same  result  could  be  obtained  by  the  order 
(10, m)  followed  by  (15,00);  however,  such  an  order  pair  not  only  requires 
twice  as  many  FLINT  instructions  but  it  also  destroys  the  contents  of  our 
regular  X  accumulator.  The  orders  (17,ni)  and  (l8,m)  invoke  both  the  X 
accumulator  and  the  Y  register,  causing  the  product  of  the  floating-point 
number  in  address  m  by  the  floating-point  number  in  the  Y  register  to  be 
added  to  or  subtracted  from  the  contents  of  the  X  accumulator,  leaving 
the  result  in  the  X  accumulator.  This  set  of  orders  permits  the  expedient 
calculation  of  the  dot  product  of  two  vectors  --  an  extremely  common  form 
of  numerical  manipulation.  -  The  (19,in)  order  is  an  inverse  division  which 
permits  the  contents  of  floating  register  m  to  be  divided  by  the  contents 
of  X.  Experience  has  shown  that  this  order  is  used  at  least  as  frequently 
as  the  normal  division  (lU,m).  It  thus  seems  very  desirable  to  have  an 
inverse  division  in  future  machines. 
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11,3  Integer  Arithmetic  Ordera 

In  most  digital  computing,  expediency  dictates  that  the  basic  program 
be  written  in  an  iterative  fashion.  That  is,  the  same  set  of  orders  is  per- 
formed over  and  over  again,  either  identically  or  with  alight  modifications 
of  the  address  parts  of  the  orders,  thereby  performing  the  same  operation 
systematically  on  large  sets  of  similar  data.  Both  in  going  around  simple 
loops  and  in  applying  the  same  set  of  orders  tc  large  sets  of  data  it  is 
necessary  to  count  how  many  times  we  have  been  around  and  to  modify  addresses. 
For  all  these  purposes,  it  -.s  necessary  to  do  simple  arithmetic  with  integers. 
In  the  IAS  machine  integers  may  be  interpreted  as  addresses  in  the  high-speed 
memory  if  they  lie  in  the  range  from  0  to  1023 •   I'^  order  to  carry  out 
integer  arithmetic,  a  second  arithmetic  organ  is  provided  by  FLINT.  The 
accumulator  of  the  integer  arithmetic  organ  is  designated  by  A  and  its 
contents,  in  accordance  with  previous  usage,  by  a.  We  have  a  set  of 
instructions  applying  to  the  A  accumulator,  which  is  quite  similar  to  our 
floating  arithmetic  orders.  We  can  Fetch  numbers  into  A,  Add,  Subtract, 
Multiply  and  Store  the  results  of  A  back  in  storage  registers.  The  detail- 
ed order  code  is  given  below. 


(00, n) 

(n) 

>  A 

(Ol.n) 

a  +  (n) 

>  A 

(02,n) 

a  -    (n) 

■>  A 

(03, n) 

a   .    (n) 

>  A 

(OU,n) 

a 

i     (n)    1          ^) 

(05,n) 

a 

■>(n) ,  non-clear  A 

We  have  again  some  thirty  storage  registers  but  they  are  not  the  same 
storage  registers  which  were  referred  to  by  the  floating-point  arithmetic 
orders.  These  storage  registers  are  used  only  for  storing  integers  in  the 
range  from  0  to  1023  (that  is,  integers  representable  by  10  binary  digits). 

The  A  accumulator  is  identical  with  the  integer  register  whose  address 
is  0.   (This  in  contradistinction  to  the  X  acciomulator,  which  is  not 
identical  with  any  of  the  addressable  floating-point  storage  registers). 
The  choice  of  integer  register  0  for  the  A  accumulator  has  been  intentional 


*■)  expla.•^^ii.(^v^  cw  r-ve,-xt  pcx^e 
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30  that  the  order  (00,00)  should  be  a  no-operatlon  order  --  since  it 
fetches  the  contents  of  integer  register  0  to  itself. 

The  operation  of  division  is  missing.   If  ve  divide  two  integers 
the  result  is,  generally  speaking,  not  an  integer,  and  hence,  this 
operation  is  omitted.  Instead,  the  (Oi^,n)  order  has  been  used  for  a 
rather  different  type  of  command,  namely  a  question  or  test  which  asks 
whether  the  current  contents  of  the  A  register  are  unequal  to  the  current 
contents  of  integer  storage  register  n?  It  is  this  question  which  is 
used  primarily  to  decide  whether  we  have  been  through  a  particular  sequence 
of  operations  a  sufficient  number  of  times  or  whether  we  ought  to  keep 
on  going  through  it  again.  This  order  is  used  only  in  connection  with 
a  "transfer"  order. 

11. Ij-  Transfer  Orders 

The  FLINT  code  is  a  sequenced  language,  that  is,  FLINT  obeys  the 
instructions  in  the  order  in  which  they  are  written  down,  unless  a  transfer 
instruction  is  encountered  which  tells  the  machine  to  go  to  some  other 
location. 

FLINT  possesses  two  transfer  orders  —  (30,w)  and(31,w).  These  orders 
behave  identically  as  far  as  their  transfer  function  is  concerned,  so  they 
will  be  discussed  together.   (The  31  orders  have  additional  properties 
which  will  be  discussed  later.)   If,  in  the  normal  sequence  of  eveiits,  our 
FLINT  code  encounters  an  order  (30,17)  it  will  seek  its  next  instruction 
at  the  beginning  of  word  17  of  the  current  line  of  orders. 

While  the  (30,w)  instruction  is  sufficient  for  transferring  to  arbit- 
rary words  within  one  line,  we  must  also  have  a  method  for  transferring 
to  arbitrary  words  within  other  lines.   This  clearly  takes  more  informa- 
tion than  we  pack  into  one  order,  and  so  FLINT  recognizes  two  kinds  of 
transfers  -  short  and  long.  For  the  long  transfer  it  is  necessary  to  specify 
to  what  line  we  wish  FLINT  to  go  next.  As  lines  of  FLINT  coding  are 
stored  on  the  magnetic  drum,  a  particular  line  of  orders  has  an  address 
consisting  of  (d)  where  d  may  be  any  integer  from  0  through  11  and 
{Jl)    =  any  integer  from  0  through  31*  A  long  transfer,  then,  consists  of  two 
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successive  Instructions  within  the  same  word,  having  the  general  form 
(30 ,w)  (d,  ,  )• 

The  need  for  long  transfers  and  the  necessity  for  transferring 
to  the  first  instruction  of  a  word  imposes  restrictions  on  the  place- 
ment of  transfer  orders  within  our  code.   In  particular,  it  requires 
that  all  transfer  orders  must  be  the  last  effective  order  in  the  word 
in  which  they  are  written.   If,  therefore,  we  write  a  transfer  instruc- 
tion in  the  first  two  positions  of  a  word,  the  last  two  positions  are 
wasted.  Fortunately,  instruction  space  is  plentiful  within  FLINT  and 
this  need  not  worry  us. 

Both  the  (30,w)  and  (31^w)  transfer  orders  are  conditional  on  the 
state  of  Transfer  Counter  (TC) .  This  counter  is  normally  in  a  YES 
condition,  but  the  (0U,n)  comparison  order  will  cause  it  to  be  set 
to  whatever  the  answer  may  be  to  the  question  it  asks.  Any  transfer 
order  will  be  obeyed  if  the  TC  is  in  a  YES  condition  at  the  time  the 
transfer  order  is  encountered  by  FLINT.   If  the  TC  is  NO,  transfer 
orders  will  not  be  obeyed  —  control  passing  to  the  first  instruction 
of  the  next  word  --  but  the  TC  is  reset  to  YES. 

11.5  Indirect  Address  Orders 

In  most  digital  computing  we  deal  with  two  rather  different  types 
of  data:  simple  isolated  floating-point  numbers  such  as  constants, 
results  of  previous  computations,  and  isolated  coefficients  --  and  with 
large  systematic  groups  of  data,  vectors  or  matrices,  all  of  which  are 
to  be  treated  in  some  systematic  fashion.  The  31  floating-point  storage 
registers,  together  with  the  arithmetic  orders  10  through  19>  suffice 
to  treat  the  isolated  individual  numbers,  but  they  clearly  cannot  handle 
large  quantities  of  systematic  data.  To  handle  these  numbers,  we  have 
the  indirect  address  orders .  They  are  very  similar  to  the  direct-address 
arithmetic  orders : 

(20,n)      ((n))  — >  X 

(51,n)  X  +  ((n))  — >X 


(29, n)  ((n))  /  X  — >X 
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except  that  ((n)),  0  <  n  <  30,  stands  for  the  floating-point  number  vhoae 
address  is  presently  stored  in  #n  of  the  "integer"  storage  location.  More 
exactly:  when  FLINT  encounters  the  order  21.17>  it  will  find  in  register 
17  an  integer  somewhere  between  0  and  1023,  which  it  will  then  interpret  as 
an  address,  and  it  will  go  to  that  address  to  find  the  floating-point 
number  which  it  will  then  add  into  the  accumulator.   In  the  midst  of  this 
process  it  will  also  transfer  the  contents  of  integer  register  17  to  A. 
These  indirect -ad dress  orders  are  at  once  the  most  useful  and  the 
most  difficult  to  use  in  the  FLINT  system,  difficult  especially  for  occa- 
sional users  who  have  not  previously  been  exposed  to  coding  for  our,  or 
similar,  machines.   If  compared  with  the  corresponding  procedure  for 
"direct  coding",  the  use  of  these  orders  is  a  simplification  in  as  far 
as 

(i)  no  distinction  of  left  or  right  hand  phases  is  necessary 
(ii)  if  the  same  address  appears  in  several  orders  and  has  to 

be  modified  in  all  of  them,  only  one  "integer"  needs  to 

be  modified  rather  than  all  order  words, 
(iii)  the  reading  of  the  address  to  be  modified  (into  A)  is  automatic 

with  the  indirect  address  order. 
Several  attempts  have  been  made  to  replace  the  indirect-address  orders 
by  another  scheme;  but  none  seemed  simpler  without  sacrificing  flexibility. 
This  seems  to  be  inherent  in  the  interpretive  system.   In  other  words, 
only  a  compiler  --  other  than  a  word  by  word  translater  --  could  bring 
real  progress.  Such  a  project,  however,  could  not  be  achieved  without 
much  more  time  and  manpower. 

11.6  Non-Addresa  Orders 

So  far,  all  orders  have  been  divided  into  an  operation  and  an  address 
part,  the  latter  one  giving  the  address  of  one  of  the  operands.  The 
remaining  orders  06.00  through  09*31  form  an  exception  to  this  rule. 
The  first  group  is  used  for  multiplying  or  dividing,  respectively,  a 
floating-point  number  by  an  integer  up  to  31,   which  is  directly  given 
in  the  second  pentad: 

(06.k)     X  .  k  >  X 

(07. k)     X  /  k  >  X 
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These  orders  have  been  useful  for  the  coding  of  formulae  containing 
small  fixed  integers,  such  as  integration,  differentiation  and  in- 
terpolation formulae,  where  they  relieve  the  coder  from  the  necessity 
of  prestoring  such  numbers  (as  2,  3>  ^>  6  etc.)  as  floating-point 
numbers . 

The  second  group  will  replace  the  number  x,  stored  in  the  X 
accumulator  by  a  function  of  x.  At  present  the  list  is  short; 

f  08.00)        "f- — >x 

foS.Ol)  loggX  >  X 

(08.02)  exp  X  >  X 

(08.03)  log^x  ^>  X 

(08. oil-)  Coshx >  X 

(08.05)  Sinh  X >  X 

but  an  extension  is  easily  possible  by  using  the  magnetic  drum  for  the 
storage  of  such  subroutines.  A  code  for  the  four  baaic  trigonometric 
functions  has  been  written  for,  but  not  yet  incorporated  in,  the  FLINT 
system,  and  an  algorithm  for  arctan  x  has  been  developed. 

The  third  group  is  a  somewhat  accidental  mixture  of  non-addresa 
orders  which  have  proven  to  be  quite  useful  for  various  purposes: 


(09.00) 

1     ^>  A 

(09.10) 

1  >  X 

(09.01) 

a  +  1  >     A 

(09.11) 

X  +  1  >  X 

(09.02) 

a  -   1  >     A 

(09.12) 

X  -  1  >  X 

1  \ 

(09.03) 

(5 .   a  >    A   ■""' 

(09.13) 

(3'.   X         >X  ^^ 

(09. oi^) 

-  a  >     A 

(09. li^) 

-  X  >  X 

(09.05) 
(09.06) 

a     <    0  >   2) 
a     =0^2^ 

(09.15) 
(09.16) 

X     <  0^^^ 
X     =0^^^ 

(09. 07) 

a     >    0  ^  '^ 

(09. IT) 

X     >  0%    ' 

'        (09.23)  will  deposit  the  signum  of  x  in  a  special  register  --  acces- 
sible only  by  (09.03)  and  (09-13)  --  and  leave  the  absolute  value 
of  X  in  X.   (09.03)  and  (09-13)  will  multiply,  without  changing  c?  ,  the 
contents  of  A  and  X,  respectively,  with  this  signum. 

^^   These  "questions"  are  used,  just  as  the  (Oi^,n)  order,  to  make  the 

next  transfer  order  --  which  need  not  be  the  next  order  —  condition- 
al. 


3) 

M 

5) 
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(09.18) 

^ 

>a3) 

(09019) 

0  < 

«>  a  +  X    1 
X     <  1  ^)    J 

(09. 2i^) 

ixi 

>  X 

(09.25) 

-1 

>  X 

(09.26) 

2 
X 

^>  X 

(09.08)  2^  .  X  ^>  X 

(09.09)  a  >  X 

a  integer 

(09.20)      0 ^>  X 

(09-21)     2  X  >  X 

(09.22)  |x  •>X 

(09.23)  .gnx  >(^  (09.27)  x^/2   ^^ 

andlx( ^>  X  -^^ 

(09.28)  reverses  the  position  of  the  "transfer  condition"  from  YES 
to  NO,  or  NO  to  YES,  respectively. 

(09.29)  is  used  only  in  connection  with  the  31  transfer  orders.  The 
only  difference  between  a  30  and  a  corresponding  31  transfer 
is  that  the  latter  will  cause  FLINT  to  "remember"  (Store) 

the  location  where  the  31  transfer  was  given.  The  next  (09.29) 
order  will  then  transfer  back  to  the  beginning  of  the  order 
word  following  the  memorized  address.  This  order,  placed  at 
the  end  of  a  subroutine,  facilitates  references  to  this  sub- 
routine at  various  places  in  the  code  (example:  evaluation 
of  f (x,y)  at  the  various  stages  of  the  RUNGE-KDTTA  procedure) 
without  the  coder  having  to  worry  about  changing  the  final 
transfer  of  the  subroutine.  This  (09.29)  transfer  is  condition- 
al in  the  same  sense  as  the  other  transfers,  if  preceded  by 
questions,  such  as  the  0U,n  or  09.O5  etc.  orders.  As  the 
subroutine  referred  to  might  itself  refer  to  a  "subsubroutine", 
care  has  been  taken  to  provide  for  up  to  four  levels :  each 
31  transfer  will  lead  one  step  "down",  each  (09.29)  one  step 
"up". 


(09.18)  brings  the  exponent  of  x  to  A  without  changing  the  contents 
of  X.  Particularly  useful  for  testing  the  order  of  magnitude  of  x. 

This  is  about  the  most  direct  transfer  possible  from  X  to  A  since 
a  must  be  integer.  Used  e.g.  to  find  the  "reduced  argument"  x 
for  a  periodic  function. 

see  footnote  (1)  on  previous  page  11.6.0 
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(09«30)  Read  in  decimal  data  from  IBM  cards.  Each  card  holds 
up  to  two  floating-point  decimal  numbers,  together 
with  their  future  addresses  and  in  a  format  suitable 
for  tabulation.   (09.30)  will  cause  conversion  of  such 
cards  placed  in  the  read  hopper;  the  floating-point 
numbers  will  go  first  to  X,  then  to  the  storage  location 
specified  in  the  "address"  columns,  while  this  address 
goes  to  A.  This  process  will  stop  if  an  empty  card  or 
an  empty  half  card  is  encountered.  FLINT  will  then 
proceed  to  the  next  order. 
(09.31)  reconverts  x  into  a  floating-point  decimal  number  which  is 
stored  in  the  form  of  a  "card  image"  together  with  the 
decimal  equivalent  of  a.  Two  such  pairs  (x,  a)  fit  on  one 
card  and  are  automatically  punched  out  as  soon  as  a  third 
(09-31)  order  is  given.  The  format  is  the  same  as  for  input 
of  numbers,  as  explained  above. 
Most  of  the  functions  carried  out  by  the  orders  (06.OO)  through  (09.29) 
could  also  be  achieved  by  suitable  combinations  of  the  other  orders;  these 
non-address  orders,  however,  add  very  much  to  the  ease  of  coding  for  FLINT. 

11.7  The  "31"  Address 

FLINT  makes  provision  for  storing  both  floating-point  and  integer 
constants  right  in  the  order  code,  interspersed  among  the  FLINT  orders. 
In  effect,  the  arithmetic  orders  invoking  a  31  address  tell  FLINT  to  look 
in  "the  next  possible  space"  --  and  use  it  as  if  it  came  from  a  regular 
register.  Thus  (01,31)  will  cause  FLINT  to  add  into  A  the  next  "order" 
as  a  (binarily  punched)  integer.   If  we  wish  to  multiply  the  contents 
of  A  by  3,  we  may  write  (03,31)  (00,03),  while  (01,31)  (02.03)  will  cause 
(A)  to  be  increased  by  2x32+3  =  67.  The  only  restriction  on  the  location 
of  the  integer  is  that  it  must  be  in  the  next  order  space  after  the  in- 
struction invoking  it. 

A  direct  floating  arithmetic  order  with  a  31  address  obviously  requires 
an  entire  word  for  its  datum.  Accordingly  FLINT  looks  at  the  next  entire 
word,  interprets  it  as  a  floating  binary  datum,  and  uses  it  as  instructed. 
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The  FLINT  order  code  then  proceeds  to  the  order  following  the  31-addreB8 
order,  with  no  skipping  if  this  next  order  is  in  the  same  word.  The  word 
which  has  been  used  as  a  datum  is  skipped  automatically,  FLINT  remembering 
that  it  is  not  a  group  of  orders .   (This  automatic  skipping  holds  true  even 
if  several  Sl-address  orders  are  used  consecutively,  requiring  several 
consecutive  constants  in  the  code.  Thus,  two  floating  arithmetic  orders 
with  31-addresses  do  not  refer  to  the  same  location,  but  to  two  successive 
locations .), 

An  indirect-address  order  with  a  31-address  will  deal  with  the  floating- 
point number  stored  at  the  address  given  in  the  next  quarter  word.  This 
means,  in  effect,  that  the  whole  Williams  memory,  as  far  as  it  is  not  used 
by  the  interpretive  code  itself,  is  "directly"  addressable  by  20-bit  orders. 
For  example,  (21,31)  (13»17)  will  add  the  floating-point  number  stored 
in  line  13,  column  17,  of  the  Williams  memory  to  X  and  (25,31)  (13,17)  would 
store  (X)  at  that  position. 

It  seemed  scmewhat  doubtful  in  the  beginning  whether  the  advantages 

which  are  afforded  by  the  address  "31"  would  be  important  enough  to  justify 

*) 
the  logical  complications  involved;  'but  after  a  year  of  use  and  experience 

this  question  can  certainly  be  answered  in  the  affirmative.  As  FLINT  allows 
for  the  storage  of  only  30  integers  (a.ddresses)  and  thirty  directly  address- 
able floating-point  numbers,  storage  place  would  soon  be  at  a  premium  with- 
out the  "31"  addresses;  while  with  this  facility  it  has  proved  sufficient 
for  much  bigger  problems  than  FLINT  was  originally  designed  for.  Another 
aspect  is  equally  important:  the  possibility  of  including  constants  right 
in  the  order  code  is  a  great  help  in  the  use  of  subroutines  by  cutting  down 
their  storage  specifications  to  a  minimum. 


The  additional  code  for  all  three  applications,  i.e.  for  fixed  point, 
floating-point,  and  indirect  address  orders  adds  up  to  only  10  order 
words . 
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11.8  Illuatratlve  Examples  of  FLINT  Coding 

(i)  Computation  of  a  Table  for  J^{x) ,   x  =  0(0.1)1  . 

Bessel  functions  for  small  arguments  are  best  computed  from  their 
respective  power  series: 


^     T=Q   ki  (p+k)  t^  '^  ' 


Five  terms  are  sufficient  for  full  accuracy  if  x  <  1. 

WORD  CODE  OPERATION  STOBAGE 

0  10.01  EX"  X  in  00.01 

(|)^  to  00.02 


10.01 

E  X 

09.22 

.1/2 

09.26 

x2 

15  "02 

S(|)^ 

07.28 

■^28 

09.12 

-1 

13-02 

■ir 

07.10 

-f-18 

09.11 

+1 

13.02 

.(|)' 

07.10 

fio 

08.12 

-1 

13.02 

.(|)' 

07.  oi^ 

4^ 

09.11 

+1 

13.02 

.(|)' 

13.01 

.X 

07.12 

712 

09.31 

EEC  J^Cx) 

10.01 

B  X 

09.31 

EEC  X 

09.18 

9">  a 

09.05 

A  >  0   ? 

06.10 

.10 
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WORD         CODE  OPERATION  STORAGE 


CODE 

Oi-KKATION 

09.11 

+1 

07.10 

fio 

15.01 

S  X 

30.00 

Tr 

(ii)  Find  a  mod  b,  a  and  b  integers 

00.01  R  a  a  in  00.01 

00.00  —  b  in  00.02 

00.00  —  a  mod  b  to  00.03 

00.00  — 

05.03  S  a  mod  b 

02.02  -b 

09.07  a  >  0  ? 

30.01  Tr  (to  beginning  of  this  word) 


2 
(iii)  Solution  of  the  Quadratic  Equation  x  +  px  +  q  =  0 


09.30 

convert 

09.22 

1/2 

09.1J+ 

.(-1) 

15.03 

S    -p/2 

09.23 

sgn  X  =>  (5" 

09-26 

2 

X 

12.01 
09.27 

1/. 

X     ' 

09.13 

^X     =>    X 

11.03 

+   (-p/2) 

09.31 

REC  x^ 

19.02 

l:q 

09.31 

REC  X 

q  in  00.01 
p  in-  00.02 
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(Iv)  Square  Root  of  a  Complay  Number 
1/2 


a  +  ib 


p  g.  c(    give^u     ^     a  ir  b    =  '2 


(P  +  iq) 
Algorithm: 


=>    c 


if  p  >  0:      c   =>  a,  q/2c   =>  b 
if  p  <  0:      c  =>  b,  q/2c   =>  a 


WOBD  CODE 

0  10.01 
09.23 
09.26 
16.03 

1  17.03 
09.27 

09.13 
11.01 

2  09.2if- 
09.22 
09.27 
15.02 

3  15. Oi^ 
09.21 
19.03 
09.00 

k  09.03 

01.31 
00.03 
20.00 


OPERATION 

E  p 

■  gn  X  =>  6" 

x2 

Q  q 

1/2 

X    ' 

,   sgn  p 

+  P 

|x| 

.    1/2 

x^/2 

S  c  as   "a" 

S  c  as   "b" 

.2 

Iq 

1  =>  a 

«   sgn  p 

/a  +  3  \ 

.=>  a     / 

S  q/2c  as    "a"   or   "b" 

STORAGE 

p  in  00.01 
q  in  00.03 
a  to  00.02 
b  to  00.04 


11.9-0 


11.9  Storage  Space  and  Operation  Times 


Every  interpretive  subroutine  is  bound  to  reduce  the  storage  space 
available  for  data  and  the  speed  of  operation.  This  is  the  price  paid 
for  the  ease  in  coding.  We  vish  to  shov,  in  this  last  paragraph,  how 
high  the  price  of  using  FLINT  is  in  these  resj^cts. 

(i)  Storage;  Let  us  again  divide  the  102^)-  words  of  the  Williama 
Memory  into  32  lines  (0  through  31)  of  32  words  each,  which  actually 
corresponds  to  the  picture  displayed  on  the  screens  of  the  memory  tubes: 
Line  31  is  reserved  for  the  "integers"  as  erplained  in  11.3; 

(31.00)  =  A; (31.31),  which  according  to  11. 7  cannot  be 
directly  addressed,  holds  the  present  order  word. 
Line  30  holds  that  line  of  pseudo  code  which  has  last  been  called 

(from  the  drum  by  a  "long  transfer")  for  execution.   (30.31) 
holds  the  negative  sum  of  (30.00)  through  (30.30)  formed 
automatically  by  the  FLINT  read-in  code.  Checking  of  the 
sum  of  the  whole  line  --  which  must  then  be  zero  —  is 
automatic  with  each  long  transfer. 
Lines  2ii--29  constitute  the  main  part  of  FLINT,  interpreting  and 

executing  all  orders  except  (O6.OO)  through  (09.31). 
Lines  21-23  are  normally  used  for  the  execution  of  the  orders  (09.OO) 
through  (09.29);  if,  however,  a  conversion,  (09.3O)  or 
(09.31),  is  ordered,  the  respective  subroutine  --  exactly  3 
lines  —  will  be  called  from  the  drum  to  these  3  lines 
of  the  Williams  Memory.  The  time  lost  by  these  drum  to 
Williams  transfers  is  negligible  compared  with  the  IBM  card 
reading  or  punching  times  which  are  necessarily  connected 
with  these  conversion  orders. 
Lines  18-21  are  normally  used  by  the  functional  subroutines  (O8.OO) 
through  (08.05);  but  if  these  are  not  needed  in  a  code,  the 
space  may  be  used  for  data. 
Lines  O-I7  are  free  for  floating-point  data,  but  only  line  0  can  be 
addressed  by  the  direct  floating-point  orders. 
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Lines  l8-2l  are  normally  used  by  the  functional  subroutines  (08.00) 
through  (08.05);  but  if  these  are  not  needed  in  a  code, 
the  space  may  be  used  for  data. 
Lines  0-1?   are  free  for  floating-point  data,  but  only  line  0  can 

be  addressed  by  the  direct  floating-point  orders. 
In  the  present  version  of  FLINT,  the  magnetic  drum  can  only  be  used 
for  the  storage  of  code,  but  not  for  data  without  referring  to  "direct 
coding"  for  the  transfer  of  such  data  to  and  from  the  drum.  There  is  no 
good  reason  for  this  restriction  and  a  corresponding  change  of  FLINT  is 
planned . 

(ii)  Operation  Times;  The  interpretation  and  execution  of  a  floating- 
point order  in  FLINT  takes  approximately  5  msec,  while  fixed  point  and  non- 
address  orders  normally  take  a  little  less  (3  to  h   msec).  For  comparison, 
machine  multiplications  and  divisions  last  approximately  one  millisecond, 
addition  and  similar  operations  a  little  less  than  0.1  msec.  It  thus  seems 
that  the  use  of  FLINT  would  lengthen  machine  times  by  a  factor  between  5 
and  50.  Better  estimates  may  be  derived  from  the  following  considerations: 

(a)  The  number  of  "mathematical"  additions,  i.e.  excluding  those  for 
bookkeeping  and  address  modification,  is,  in  general,  not  too 
different  from  that  of  the  multiplications  and  divisions  together. 

(b)  The  number  of  orders  required  to  carry  out  a  typical  sequence 

of  operations  is  normally  much  less  in  FLINT  than  in  direct  coding, 
even  if  fixed  point  coding  is  feasible,  the  average  being  about 
1:2  (excluding  FLINT  subroutines  such  as  08.OO,  09*27,  09-30, 

09.31). 

(c)  The  greater  flexibility  of  the  FLINT  order  code,  especially  the 
absence  of  scaling  difficulties,  often  allows  the  use  of  more 
efficient  mathematical  methods,  especially  since  coding  spac* 
is  practically  unrestricted. 

(d)  The  speed  of  in-and  output  is  basically  the  same  in  FLINT  and 
in  direct  coding;  hence,  if  the  time  required  for  card  reading 
and  punching  is  an  appreciable  fraction  of  the  total  running 
time  of  the  code,  the  time  factor  is  further  reduced. 
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(e)  Debugging  time«  for  FLIKT  codei  are,  even  with  the  inefficient 
present  "TRACER",  normally  much  shorter  than  in  direct  coding. 
Fewer  errors  are  made  and  those  made  are  more  readily  detected. 
The  result  can  be  summarized  as  follows.  The  average  effective  slow- 
down factors  are  approximately: 

Considering  (a)  alone 1^ 

(a)  and  (b) 10 

"     (a)  through  (c) 5-10 

"     (a)  through  (d) 2-10 

"     (a)  through  (e) 0.1-10 

The  last  figure  shows  that  the  use  of  a  floating-point  interpretive 
routine  need  not  always  be  wasteful  of  machine  time.  Experience  has  con- 
firmed that  it  can  even  lead  to  a  substantial  saving. 


12.0 

12.        ANOTHER  FLOATIHG-POIMT   lUTERFRETrTE  ROUTIWE    (Code   I906) 
Originator  and  Coder;      Irving  N.  Rabinowitz 

12.0  Introduction 

1906  is  a  floating-point  interpretive  routine,  similar  in  some  res- 
pects to  FLINT,  but  basically  different  in  its  logical  structure.  However, 
some  parts  of  the  two  interpreters,  notably  the  arithmetic  operations  are 
identical,  and  so  will  not  be  discussed  in  great  detail. 

1906  will  accept  a  pseudo-code  residing  in  the  machine,  and  will  inter- 
pret this,  executing  both  arithmetic  and  logical  operations.   In  addition, 
the  structure  of  1906  allows  the  coder  to  mingle  direct  code  with  his  inter- 
preted code  with  a  minimum  of  effort.  However,  the  orders  available  to  I906 
are  extensive  enough  so  that  this  need  not  be  done  except  to  execute  drum  read 
or  write  orders,  which  are  not  provided  for  in  the  pseudo-code. 

The  pseudo-code  for  I906  is  kept  in  the  Williams  memory  at  all  times, 
and  thus  can  be  modified  by  itself.   In  this  respect  the  pseudo-code  is  closer 
to  machine  code  than  in  the  FLINT  routine,  the  difference  being  that  between 
a  sequence -controlled  machine  and  a  stored-program  machine. 

12.1  Structure  of  the  Pseudo-machine 

As  mentioned,  I906  is  a  stored  program  machine  (and  we  might  as  well 
adopt  the  view  that  it  is  a  pseudo-machine,  since  it  avoids  circumlocutions). 
The  programmer  has  available  to  him  without  restriction  672  words  (=21  lines) 
of  high-speed  storage,  which  can  hold  numbers,  orders,  or  pseudo-orders. 
Numbers  may  be  of  two  kinds,  floating-point  numbers,  which  represent  quantities 
actually  used  in  the  computation,  and  fixed-point  numbers,  or  integers,  which 
are  used  for  logical  manipulations. 

Since  the  arithmetic  portion  of  the  pseudo-machine  is  identical  to  that 
of  FLINT,  we  need  not  describe  the  structure  of  floating-point  numbers  (see 
section  11.1) . 

Fixed  point  numbers,  or  integers,  may  be  real  numbers,  i.e.  signed  forty 
bit  machine  numbers,  which  may  be  manipulated  by  a  set  of  pseudo-orders,  but 
which  are  better  handled  by  direct  coding;  or  they  may  be  lO-bit  numbers  with- 
out sign  for  specifying  certain  logical  parameters  such  as  indices. 

The  structure  of  the  pseudo-order  words  is  similar  to  that  of  machine 
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oTder  words  in  that  two  pseudo-orders  occupy  a  single  word,  but  the  structure 
of  an  order  is  not  (address,  order),  but  rather  (address,  tag,  order).  The 
address  is  a  ten-bit  number  representing  addresses  (00,00)  through  (20,31) • 
The  tag  is  a  five-bit  number  whose  individual  bits  specify  certain  B-registers. 
The  order  is  a  five-bit  number.  This  allows  a  total  of  32  orders,  but  by 
sacrificing  the  property  of  "B-modif lability"  for  some  orders,  namely  those 
with  order  part  31,  it  is  possible  to  increase  the  number  of  available  orders 
to  63,  which  is  sufficient  for  all  purposes. 

12.2  Arithmetic  Orders 

The  arithmetic  unit  of  1906  consists  of  the  pseudo-registers  R  and  Q 
both  of  which  hold  floating-point  numbers.  The  orders  affecting  these  reg- 
isters are  the  following:   (Let  (x)  mean  "the  contents  of  x",  where  x  may 
be  a  memory  location, or  E  or  Q. 


A) 


B-modif iable: 

2h 

Bring 

05 

Add 

15 

Subtract 

25 

Subtract  Absolute  Value 

06 

Multiply 

16 

Multiply  negative 

07 

Multiply  8e  accumulate 

17 

Multiply  &  accumulate 
negative 

08 

Divide  E 

18 

Divide  x 

26 

Store 

27 

Store  zero 

29 

Load  Q 

(x)  — >  R 
(E)  +  (x)  — >  E 
(E)  -  (x)  — >  E 
(E)  -|(x)|— >E 
(E)  -f  (x)  — >  E 
-(E)  X  (x)  — >  E 
(E)+(Q)(x)— >R 

(E)-(Q)(x)— >E 

(E)  /   (x)— >  E 
(x)  /   (E)— >E 
(E)       — >  X 
0  — >  x;  E  unchanged 
(x)— >  Q 


B)  Not  B-modlf iable,  and  address  not  used: 

0031  Square  root  nJIeY  — >  E 

0131  Exponential  e^  '   — >  R 

0331  Absolute  value  |(R)|  — >  E 

OU31  Negative  -(E)  — >  E 

0631  Square  (E)^  — >  E 
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C)  Not  B -modifiable,  address  used; 

0531  Multiply  (E)  by  address  (E)  xN  — >  E 

0731  Add  (E)  and  address  (E)  +  N  — 5€? 

(In  these  two  orders,  N  is  one-fourth  the  address,  so  that  the  numbers 
from  l/k   to  127  3A  are  available  for  multiplication  and  addition  without 
being  stored.  Negative  numbers  may  also  be  used  by  complementing  the  address). 

The  nature  of  these  arithmetic  orders  is  self-explanatory,  so  we  need  not 
go  into  them  further,  except  to  note  that  they  exist  in  the  two  classes  of  B- 
modifiable  and  not  B-modifiable.  These  terms  will  be  explained  when  we  describe 
the  logic  of  1906. 

12.3  Fixed-point  Arithmetic 

In  addition  to  the  pseudo-registers  E  and  Q,  there  exists  another  pseudo- 
register  called  C,  in  which  fixed-point  numbers  may  be  m.anipulated.  The  orders 
affecting  C  are  not  B-modifiable,  but  do  use  the  address  portion  to  represent 
a  storage  location.  The  orders  are 

0831  Bring  (x)  — >  C 

0931  Bring  negative  -(x)  — ■>  C 

1031  Add  (C)  +  (x)  — >  C 

1131  Subtract  (C)  -  (x)  — >  C 

1231  Multiply  2"^^(C)(x)  — >  C 

1331  Store  (C)  — >  X 

These  are  all  fixed -point  operations  except  for  1231,  Multiply,  which^ 
as  in  FLINT,  includes  a  left  shift  of  9  places.  Since  manipulations  on  fixed- 
point  numbers  are  primarily  for  logical,  rather  than  arithmetic,  purposes, 

9 
the  fixed-point  operations  will  be  applied  to  integers  modulo  2  ,  which  in 

a  machine  word,  of  course,  are  represented  as  fractions*  So  that  a  product 

be  properly  scaled,  the  left  shift  of  9  places  is  included  after  multiplication. 

12. i|-  Transfer  or  Jump  Instructions 

In  order  to  change  the  sequence  of  operations,  several  types  of  trans- 
fers have  been  included  in  the  order  code.   It  is  possible  to  change  the 
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sequence  of  operations  unconditionally  by  the  order 

10  Jump  Take  next  order  from  address  and  phase 

specified. 

Note  that  jumps  may  be  made  to  either  order  in  a  word.  This  specification 
is  accomplished  by  a  phase  bit,  which  is  part  of  the  tag.   In  jump  orders  this 
bit  is  interpreted  as  specifying  the  feft-hand  order  of  the  location  given 
by  the  address  if  the  bit  is  zero,  or  the  right-hand  order  if  it  is  one. 

The  program  may  be  stopped  by 

00  Stop  and  Jump  Stop^   If  restarted,  execute  10.  However 

if  the  address  is  (00,00),  go  on  to  the 
next  order. 

Certain  jumps  may  be  conditional  on  the  state  of  registers.  Two  such  are  con- 
ditional on  the  state  of  R,  two  on  the  state  of  C,  and  l6  on  the  state  of 
l6"breakpoint  switches". 

The  orders  conditional  on  the  state  of  E  are 

20  Positive  conditional  jump   If  (E)  ^  0,  execute  10.   If  (E)  <  0, 

no  operation 

30  Negative  conditional  jump   If  (B)  <  0,  execute  10.  Otherwise 

no  operation. 

It  may  be  noted  that  there  is  no  jump  conditional  on  the  state  (E)  =  0. 
The  reason  for  this  is  that  in  floating-point  arithmetic,  round-off  effects 
can  cause  troubles  in  this  respect. 

The  two  orders  conditional  on  the  state  of  C  are; 

li<-31  Jump  if  (C)  is  positive  or  zero. 
1531  Jump  if  (C)  =0. 

Note  that  here,  since  we  are  dealing  with  integers,  a  zero  test  is 
included,  and  is  indeed  a  very  powerful  method  of  terminating  loops,  if  the 
programmer  is  counting  his  way  through  a  loop,  a  process  which  in  general 
is  unnecessary,  due  to  the  power  of  the  B-registers. 

The  other  sixteen  conditional  jump  orders  are  variants  of  the  order 

23  Breakpoint  stop  and  jump.   If  "breakpoint  switch"  n  is  1,  execute 

00;  if  it  is  0,  no  operation. 

Ideally,  of  course,  the  breakpoint  switch  would  be  a  true  mechanical  switch 
which  the  operator  could  turn  to  an  "on"  or  "off"  position.  However,  what 
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is  actually  done  is  to  reserve  a  word  in  storage  for  this  purpose.  The 
first  sixteen  bits  of  this  word  are  the  l6  breakpoint  switches,  and  are 
manually  made  1  or  0  by  the  oi^erator,  at  the  instruction  of  the  programmer. 
When  the  interpreter  encounters  the  order  (a,0,n,23)  .■>  it  inspects  bit  n  of 
the  breakpoint  word.   If  it  is  a  one,  the  order  00  is  executed,  i.e.   it 
stops,  and  upon  restarting  executes  a  jump  to  the  order  at  phase  ^  in  word 
a.   If  it  is  a  zero,  the  interpreter  ignores  this  order  aad  proceeds  to  the 
next  order.  The  purpose  of  this  order  is  to  aid  in  debugging  a  program.  The 
programmer  can  insert  23  ordera  at  various  poiats  '7f  hic(  prog^^am^ which  will 
transfer  to  a  prlat-out  routine  wiiicli  he  feass  "sirit-ieR.  Wikfea  fii's-^  ruaning  his 
program,  he  can  set  all  tiie  'b:-eakyoint  bits  xo  onts,  so  t'jMx   -she  pro^'aa  will 
stop  at  each  23  order.  Af'cer  the   pr:'.?a--out ,  --he  fiv'P!-';  'jf:',   in  tliis  breakpoint 
word  is  made  ssro,  and  the  progi'am  oontinued.  ¥xh  ■'!;h'i  Tm.nhi.vie  will  not  stop 
at  the  order  00P3,  but  will  stop  at  the  ne^ri  'hi'e&h'^joia-c-  stop  and  ju.ap  encounter- 
ed. The  process  can  then  be  repeated  with  all  the  breakpoiats  in  the  program, 
thus  gii'ing  the  progranmer  lnforri»e,tion  as  t.i  th*?  jtat'?  o^'  the  memory  at  inter- 
mediate points  of  the  computati'-in.  When  the  progv^am  is  finally  debugged,  it 
need  not  be  changed  at  all,  sin-:e  the  23  orders  ha^e  ao   effect  on  th^  execution 
of  the  correct  progj:'am,  the  breakpoint  swi-<iches  not  be^.-ag  set. 

12 . 5  Other  Logical  Instructions 

T^To   other  logical  ius-JTuctioua  are  aveiiJ.abie  to  the  progi^ammer s 

21  Set  return  If  this  order  is  encountered  in  word  a, 

then  in  the  phase  and  address  specified, 
there  is  written  the  half-word  <  a+1,  10  > 

22  Exit  Exit  from  interpretive  mode  to  direct 

coding. 

Order  21  is  designed  for  use  with  closed  subroutines.  Symbolically, 
its  use  is  as  follows :  when  a  closed  subroutine  is  to  be  used  from  word  a, 
the  programmer  writes  <Exit,  2l/Entry,  10  >,   the  effect  being  that  when  the 
closed  subroutine  has  been  executed,  and  control  has  arrived  at  the  exit  word 
of  the  subroutine,  a  jump  is  made  back  to  the  main  program.. 

The  order  22  is  used  to  go  out  of  the  interpretive  mode  into  direct  cod- 
ing, and  has  the  form  of  a  10  order,  except  that  the  order  at  the  specified 
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address  and  phase  is  not  interpreted,  but  is  executed  by  the  real  machine. 
To  re-enter  the  interpreter,  the  usual  tjrpe  of  link-planting  order  pair  is 
used,  i.e.,  to  re-enter  from  word  a,  the  programmer  writes  <  a  CP/3100  Ul  >, 
in  direct  machine  code. 

12.6  The  B-registers  and  their  Instructions 

We  now  come  to  the  most  powerful  part  of  code  I906,  the  B-registers. 
These  are  a  set  of  registers  labelled  Bl,  B2,  B^,  and  b8  (or  Bl,  B2,  B3,  bU, 
if  desired.  The  1,  2,  k,   8  notation  is  a  more  mnemonic  one,  as  will  be  aeeaj , 
which  have  the  capabilities  of  counting  and  modifying  addresses,  30  that  loop 
formation  is  reduced  to  the  writing  of  two  instructions.  It  will  be  noted  that 
the  arithmetic  orders  (  §12.2  )  were  divided  into  several  groups,  one  of  which 
is  that  of  the  B~modifiable  instructions.  For  these  instructions,  the  tag  part 
of  the  order  playa  a  role.  The  tag  consists  of  five  bits,  four  of  which  refer 
to  B-registers.   It  has  the  form  0  bn  b,  b  b  ,  where  b.  refers  to  Bi.  The 
reference  is  the  following:  if  b.  is  one,  the  address  used  by  the  pseudo- 
machine is  not  the  address  written  in  the  order,  but  is  this  address  plus  the 
contents  of  Bi.  This  is  the  effective  address.  The  use  of  this  device  is  to 
allow  references  to  arrays  of  numbers  without  having  to  go  through  the  process 
of  counting  and  address  modification  within  a  loop.   If  more  than  one  b  is  one, 
the  effective  address  is  the  original  (written)  address  plus  the  contents  of  all 
the  Bi  referred  to  by  the  b  .   In  this  way  it  is  extremely  convenient  to  refer 
to  two-  or  three-dimensional  arrays. 

Besides  their  address -modification  properties,  the  B-registers  have  count- 
ing properties  which  are  used  to  go  around  loops  the  proper  number  of  times,  and 
to  exit  from  the  loop  when  it  has  been  executed.  This  is  accomplished  by  the 
use  of  two  types  of  orders,  the  SET  and  TEST  orders.  These  are 

01  Set  Bl  Each  of  these  is  a  three-address  order,  affecting 

the  B-register  mentioned.  These  three  addresses 

02  Set  B2  are  those  of  I,  A.  ,  and  F,  respectively,  and  have 

the  effect  of  setting  I, A,  and  F  into  three 

03  Set  B3  storage  locations,  the  totality  of  which  consti- 

tute the  pseudo-register  Bi. 
0^  Set  Bh  1   is  the  initial  value  of  an  index,  £\  is  its 

increment,  and  F  its  final  value. 


11 

Teat  Bl 

12 

Test  B2 

13 

Test  B3 

lU 

Test  Bk 
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These  orders  test  the  B-register  referred  to  in 
order  to  see  whether  the  loop  has  been  traversed 
I  (F  -  I)  /Altimes.   If  this  is  the  case,  then  con- 
trol passes  to  the  next  order  in  sequence.   If  not 
then  the  contents  of  the  B-register  are  increased 
hy  A  (which  may  be  positive  or  negative),  and 
a  jump  is  made  to  the  address  and  phase  specified, 
in  order  to  traverse  the  loop  again. 

The  use  of  these  orders  is  as  follows:  as  an  example,  let  us  consider  the 
simple  problem  of  adding  up  the  numbers  located  at  storage  locations  a,  a+1,  a+2, 
...,  a+99.  The  simple-minded  method  is,  of  course,  to  write  down  100  orders  of 
the  form  a  Bring 

a+1  Add 

a+2  Add 


a+98  Add 
a+99  Add. 

There  are  obvious  disadvantages  to  such  a  course.  The  mcare  usual  procedure  is  to 
set  up  a  loop  to  accomplish  the  same  thing.   In  machines  without  address -modify- 
ing capabilities,  this  loop  is  accomplished  by  setting  a  counter  to  zero,  using 
it  to  change  the  address  of  an  order,  executing  this  order,  testing  the  counter 
against  its  final  value,  and  either  increasing  it  and  going  back  to  the  beginning, 
or  leaving  the  loop.  In  the  liitfci'pretiv>;  mode,  thia  is  accoiaplisa6d  more  aimply 
by  the  use  of  the  B-registers.  For  the  problem  mentioned,  we  could  proceed  as 
follows;  assuming  that  B,  the  pseudo-accumulator  is  clear,  we  have 

X  SET  Bl    <D>,  <1>,  <99>         (Notation:  The  symbol  <5r>  means 

"the  location  of  x.") 

(x+l)L  ADD      <^o-^'  ^1 

(x+l)R  TEST  Bl   x+l,L. 

This  two-word  loop  does  the  following:  in  word  x,  the  SET  order  takes  the  cont- 
ents of  the  three  locations  specified,  namely  0,  1,  and  99,  as  integers  scaled 
2"^,  and  deposits  them  in  the  three  storage  locations  which  make  up  Bl,  the 
location  to  which  1=0  goes  being  used  internally  to  I906  as  the  counter  for 
the  loop.  Thus  the  first  value  for  the  index  is  taken  as  zero.  The  left  half 
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of  word  x+1  contains  the  order  to  add  the  next  value  of  a.  into  the  ac- 
cumulator. The  address  written  is  that  of  the  first  value  of  a^,  namely 
a  ,  but  the  effective  address  is  a  +  (Bl) ,  since  the  b  bit  is  in  the  tag 
part  of  the  instruction.  The  first  time  through  the  loop,  (Bl)  =  i  =  0, 
so  that  a  is  brought  to  R.  At  word  x+1,  right  phase,  is  the  order  to 
terminate  the  loop,  TEST  Bl.  The  sequence  of  events  that  occurs  here  is 
as  follows;  the  loop  counter  is  compared  with  99»   If  they  are  unequal, 
it  is  increased  by  l(the  preset  value  of  ^  ) ,  and  a  jump  is  made  to  the 
order  at  location  x+1,  left  phase,  namely,  the  ADD  order.  However,  if 
they  are  equal,  the  incrementing  and  jumping  do  not  take  place,  and  the 
next  order  to  be  executed  is  the  order  at  location  x+2,  left  phase. 

Thus  the  effect  of  the  loop  is  to  add  a  ,  a^,  .  .   .,  a   into  the 
pseudo-accumulator,  and  then  exit  from  the  loop,  with  the  sum  in  R. 

The  structure  of  such  loops  has  an  almost  exact  counterpart  in  the 
mathematical  notation  used  for  these  operations.  For  the  problem  mentioned, 
one  might  write 


..  I 


22_ 


i  =  0 


to  describe  the  formation  of  the  sum.  Another  way  would  be  to  say 

(E)  +  a.  >  R   for   i  =  0(1)99, 

which  has  its  exact  counterpart  in  the  code.  The  code,  in  symbolic  terms, 
could  have  been  written 

X    SET  i  FOR    0(1)99 
(x+l)L   ADD  a. 

(x+l)R   TEST  x+l,L 

More  than  one  B-register  may  be  used  in  a  single  loop.  Suppose  that  the 
problem  involved  is  to  add  two  vectors  a.  and  b  .  Suppose  that  the  vector 
a,  is  arranged  in  the  memory  in  locations  a^,  ^q+^>   ^o*^'    '    '    °    '  ^o'^'^'   ^'^' 
in  consecutive  locations  specified  by  <a  >  =  <a^>  +  i,  and  that  the  vector 
b  is  stored  differently,  e.g.,  in  every  other  location,  starting  at  <h^ 
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and  auch  that  (b.>  =  <b  >  +  2i.   It  ia  neceasary  to  use  two  B-regiatera 
for  address  modification  in  tlie  same  loop.  Let  us  write  the  mathematical 
discription  of  the  problem: 

c.  =  (a.  +  b.) ,  for  i  =  0(l)n 

Assuming  that  the  c.  are  stored  sequentially  starting  at  <c  >,   we  rewrite 
this  in  the  following  form: 

c.  =  a.  +  b.,  for  i  =  0(1)  n  and  i'  =  0(2)2n. 
Ill'  ^    '  ^   ' 

The  code  is  then  a  transcription  of  this  statement  to  a  vertical  format: 

X  SET  i  =  0(l)n 
x+1  SET  i'  =  0(2)2n 

x+2,L  BRING  a. 

'  1 

x+2,E  ADD  b. 

x+3,L  STCEE  c. 

x+3,B  TEST  i,  x+4,L 

X4J+,L  TEST  i',  x,2,L 

Note  that  although  there  is  only  a  single  loop,  yet  there  are  two  SET  orders 
and  two  TEST  orders,  since  two  B-registers  are  involved.  Since  the  indices  are 
to  be  counted  synchronously,  two  TESTS  must  be  made  for  each  traversal  of  the  loop. 
The  first  TEST  increases  the  index  i,  while  the  second  TEST  is  used  to  escape  from 
the  loop  when  finished »  Thus  the  address  of  the  first  is  x+i<-,L,  so  that  regardless 
of  whether  or  not  the  loop  is  done,  the  second  TEST  is  executed  whenever  the  first 
is. 

In  dealing  with  two-dimentional  arrays  such  as  matrices,  it  is  desirable  to 
be  able  to  modify  an  address  by  two  indices.  As  mentioned  above,  this  is  possi- 
ble by  including  more  than  a  single  1  in  the  tag.  Consider  the  problem  of 

the  multiplication  of  two  square  matrices.  Mathematically  we  write 
n-1 
c.^  =  ^2    a.j^b^j    for  i  =  0(l)n-l     j  =  0(l)n-I 

k=0 
However,  we  must  take  into  account  the  fact  that  the  indices  are  not  all  in- 
creased by  unity,  since  the  matrices  may  be  stored  in  the  memory  in  such  a 
way  as  to  require  different  modification.   If  we  assume  that  the  matrices  are 
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stored  in  the  memory  so  that  the  elements  of  the  matrix,  read  from  left  to 
right  and  top  to  bottom,  occupy  consecutive  locations,  e.g.,  ^qq)  ^q-\> 

'    '    '    '  ^cm-V   ^lo'  '        ''   ^l,n-l'  •  •  •  '  Vl,o'  Vl,l  '"'   Vl,n-1' 
then  a  simple  description  of  the  location  of  any  element  is 

<^ij>  =  <^oo>  -^  ^  -^  ^J 
Thus,  for  the  three  matrices  involved,  we  have 


<a 


il^ 


=  <a  >  +  i  +  nk 


=  <a  >  +  i  +  k' 
oo 


<^kj>  =  <^oo>-^^^  J* 

<*'ij>  =  <%o>  -^  ^  -^  J'  ' 

and  thus  there  are  four  indices  which  must  be  used,  namely,  i,  j',  k,  and  k', 

2  2 

such  that  i  =  0(l)n-l,j'  =  0(n)n  ,k  =  0(l)n-l,  and  k'  =  0(n)n  .  We  may  then 

write  the  code  as 

X  SET  i  =  0(l)n-l 
x+1  SET  j'  =  0(n)n^ 
x+2,L  BRING  zero 
x+2,E  BRING  zero 


(This  order  has  the  effect  of  clearing 
R  prior  to  accumulating  products.  The 
second  BRING  is  just  a  skip  in  order 
to  get  into  the  proper  phase.) 


x+3    SET  k  =  0(l)n-l 
x+U    SET  k'  =  0(n)n^ 


x+5,L  LOAD  L       a^j^, 

x+5,R  MULTSeACC       b^ 

x+6,L  TEST     k,   X46,R 

x-f6  ,E  TEST     k ' ,  x+5  ,L 

x+7,L  STORE         c^,, 

x+7,R  TEST     y,      x+2,R 

x+8,L  TEST     i,      x+l,L 
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Tliis  code  does  the  following;  the  orders  at  x  and  x+1  set  up  the 

B-regiaters  for  i  and  j'  for  the  i  and  j'  loops,  and  the  orders  at  x+7,R 

and  x+8,L  close  these  loops.  The  effect  is  that  the  inner  loop  (j')  runs 

across  the  row  for  which  the  index  i  has  a  certain  value o  When  the  end  of 

the  row  is  reached,  as  signalized  by  a  failure  of  the  TEST  at  x+7,E,  the  i 

loop  is  tested.   If  the  loop  is  not  done,  the  row  index,  i,  is  increased 

and  the  next  row  of  c.   is  computed.  Wichin  these  two  outer  loops  is  a 

single  inner  loop,  which  however,  involves  the  use  of  two  B-registers,  one 

to  count  k  and  the  other  to  count  k',  in  a  fashion  similar  to  the  problem 

of  vector  addition  given  above.  Prior  to  entering  th-i  loop,  E  is  set  to 

zero  by  the  orders  at  x+2.  The  reason  that  two  auc'a  orders  are  given  is 

that  the  SET  order  occupies  a  full  word  so  that  we  must  skip  an  order  after 

the  first  BEING  to  get  to  the  left  half  of  the  next  word.  The  inner  loop 

consists  of  the  setting  for  both  k  and  k',  the  actual  computation  of  the 

(i,j)  element  of  the  result  matrix,  the  closing  of  the  k  and  k'  loop,  and  the 

storing  of  c.  .  ^.The  actual  computation  of  c .  .  takes  place  at  x+5,  where  the 
ij  _K —  IJ 

partial  sum    >     a  ,  b,  ,  is  accumulated,  i.e.,  the  vector  product  of  row 
A— ^    Ik  Kj 

i  of  a  by  column  j  of  b.  Note  that  hsre  again,  the  TEST  at  x+6,L  is  used 

only  for  stepping  the  k  index  up,  while  the  actus-l  testing  for  end-of-loop 

is  done  at  x+6,E. 

12.7  Integer  Manipulations  and  Non-Automatic  Modification 

It  may  be  noted  that  in  the  examples  in  §12.6  we  were  dealing  only  with 
indices,  while  addresses  were  given  in  the  word  referring  to  the  operand.   It 
is  entirely  possible  to  use  the  B-registers  for  carrying  actual  addresses,  so 
that  the  code  itself  need  have  no  addresses  of  operands  in  it  at  all.  This 
is  sometimes  a  convenience,  for  example,  in  changing  codes  to  accomodate  larger 
arrays  than  was  originally  contemplated.  Thus,  if  a  programmer  has  been  using 
5x5  matrices  stored  sequentially,  and  has  to  go  to  67:6   matrices,  the  code  it- 
self would  have  to  be  changed  if  the  addresses  written  in  it  were  the  locations 
of  the  first  elements  of  the  5x5  arrays.  This  difficulty  may  be  avoided  by 
the  following  subterfuges  rather  than  dealing  with  indices  and  considering 
expressions  of  the  form  i=  l(A)  F    we  can  deal  with  addresses  in  expressions 
of  the  form  i  =  A+1(Z\)A+F,  where  A  is  an  address  which  is  stored  in  some 
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location.  The  loop  setup  would  then  involve  forming  the  integer  sums  A+I 
and  A+F  where  A,  I,  and  F  are  preset  integers^   and  SETting  the  B-registera 
with  these  numbers,  rather  than  with  the  values  of  the  indices  alone.  The 
reference  to  the  operand  a.  would  then  be  (say)  ADD  0,  b^  rather  than  ADD 
<a  >,b,,  and  it  would  become  unnecessary  to  change  the  code  itself  to  make 

0    1 

changes  in  the  placing  of  operands. 

Occasionally  it  may  be  advantageous  to  modify  addresses  in  a  coded 
fashion,  rather  than  using  the  automatic  features  of  the  B-registers,  for 
example,  if  one  were  dealing  with  expressions  which  might  involve  the  use 
of  more  than  four  B-regisi:ers  at  once.  Such  might  be  the  case  in  computing 
an  expression  of  the  form 


for  i,  j  =  0(l)n,  and  m  =  0(1)M.  All  four  B-registera  would  be  tied  up  dur- 
ing the  inner  loops  of  the  computation,  ao  that  it  would  be  necessary  to  have 
a  separate  counter  for  m,  which  could  be  added  into  the  settings  of  the  B-reg- 
isters during  the  inner  loops,  i.e.  one  could  compute  the  numbers  A+I+(M+l)m 
and  A+F+(M+1)M  before  setting  the  B-regiaters  in  the  inner  loops.  Alternatively, 
the  actual  orders  could  he  modified  by  adding  (M+l)m  to  their  addresses  prior 
to  setting  the  B-registers.  This  can  also  be  accomplished  by  means  of  the  fixed- 
point  instructions^  and  the  jumps  o  rnditionai  on  the  state  of  C  would  be  used 
to  determine  when  the  m-loop  has  been  completed. 

12.8  Decimal  Inpt\t  and  Output 

There  are  two  orders  which  affect  the  input  and  output  of  floating-point 

number a , 

09   Eead  Read  n  cards,  each  with  two  words,  into  locations 

specified  by  address,  sequentially,  n  given 
by  the  tag,  1<  n  <  32. 

19   Punch  Same,  in  the  other  direction. 

The  card  format  consists  of  two  words  per  card,  each  half  of  the  card 
containing  a  single  floating-point  number.  The  format  of  a  number  is  almost 
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completely  free,  with  the  exception  that  the  sign  of  the  number  must  occupy 
the  first  column  of  the  field,  and  all  numbers  must  have  a  decimal  point. 
Thus,  for  example,  to  input  the  number  3.1^^159,  it  is  not  necessary  to  punch 
it  as  31^1590000  01,  but  it  may  be  punched  directly  as  3.l4l59,  providing 
only  that  the  first  column  of  the  field  is  blank,  for  the  sign.  To  input 
both  mantissa  and  exponent,  e.g.,  the  number  3x10 ■'■°,  it  may  be  punched  as 
5^3. XIO,  where  <fi   indicates  a  blank  column.  This  notation  indicates  the  mantissa 
as  3,  terminated  by  a  decimal  point,  and  the  XlO  represents  the  exponent. 

On  output  the  number  appears  in  a  standard  fashion,  namely  in  the  form 
+  3-1^1590000X00.  Output  numbers  may,  of  course,  be  later  re-used  as  input, 
since  the  output  format  is  acceptable  to  the  input  routine. 

12.9  Tracing  Mode 

In  order  to  help  debug  interpretive  programs,  a  tracing  mode  has  been 
added  to  I906,  and  is  operated  by  using  a  special  card  preceding  the  deck. 
This  card  contains  two  addresses  punched  by  the  user,  the  address  at  which 
tracing  is  to  start,  and  the  address  at  which  it  is  to  end.  The  card  also 
includes  certain  standard  coding  which  puts  1906  into  the  tracing  mode  (by 
changing  a  transfer  instruction  within  I906) .  The  pseudo-code  is  executed 
at  full  speed  (except  for  a  few  orders  involved  in  the  testing  of  addresses) 
until  the  starting  address  for  tracing  is  encountered,  at  which  point  each 
instruction  executed  is  accompanied  by  a  card  specifying  the  contents  of  R 
and  of  the  B-register  involved,  if  any,  and  also,  of  course,  the  location 
of  the  order  itself.  This  procedure  is  carried  on  until  the  location  counter 
of  the  interpreter  becomes  larger  than  the  final  address  specified,  at  which 
point  tracing  is  stopped  and  the  code  continues  to  be  executed  at  full  speed. 

In  this  way,  it  is  possible  to  debug  sections  of  a  program  by  tracing 
the  earlier  portions  of  a  program  first,  verifying  that  they  are  correct,  and 
then  tracing  later  portions  without  being  forced  to  repeat  the  trace  of  the 
earlier  parts,  which  saves  a  great  deal  of  machine  time. 
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13.0  SERVICE  BOUTIHES  FOR  DECIMAL  HJ-  AND  OUTPUT 


In  order  to  help  save  time  in  the  coding  and  debugging  of  problems, 
several  routines  of  general  utility  were  written.  These  routines  were 
designed  to  be  used  in  conjuncticn  with  the  programmer's  own  code  to  help 
out  in  the  processes  of  input  and  output  of  decimal  numbers,  and  in  the 
detection  of  errors  in  a  code.  The  routines  described  here  were  written  at 
different  times,  and  are  therefore  quite  independent  of  each  other.   It  is 
planned  to  rewrite  them  in  an  integrated  form  so  that  any  service  routine 
can  be  called  in  a  standard  manner  with  a  minimum  of  special  handling. 

13.1  GEMEBAL  DECIMAL  INPUT  ROUTmS 

Originator  and  Coder:  Irving  N-  Rabinowitz 

This  is  a  short  routine  which  allows  a  programmer  to  dispense  with  most 
of  the  "bookkeeping"  involved  in  decimal-binary  conversion.   It  is  written  in 
the  form  of  a  closed  subroutine  which  uses  three  pre-set  parameters  for  the 
determination  of  the  arrangement  of  digits  in  the  cards.  Each  time  the  sub- 
routine is  entered,  it  has  the  effect  of  reading  the  next  number  from  the 
card,  translating  it  into  binary,  and  delivering  it  to  the  main  program  in  the 
accumulator.  The  three  parameters  necessary  for  the  routine  are:  D,  a  "stencil" 
or  "mask"  word  which  has  one's  in  those  stages  where  digits  are  punched  on  the 
card,  and  zero's  elsewhere;  S,  a  similar  word  which  specifies  the  columns  where 
signs  appear  on  the  card;  and  C,  the  number  of  cards  of  input.  As  an  example, 
if  the  card  format  consisted  of  five  seven-digit  numbers  and  signs,  filling  all 
forty  available  card  columns,  we  would  have 

D  =  0111111101111111011111110111111101111111, 

and 

S  =  1000000010000000100000001000000010000000. 

The  fields  of  each  card  need  not  be  the  same  size,  but  all  cards  must  have  the 
same  format.  While  the  routine  does  not  count,  it  does  signal  the  end  of  the 
input  by  the  state  of  the  Q-register.   If  this  is  non-zero,  then  the  input  is 
finished,  otherwise  there  is  more  to  come. 

Since  it  is  written  in  the  form  of  a  closed  subroutine,  the  programmer 
has  complete  freedom  as  to  the  placing  of  the  input  in  the  memory,  the  scaling 
of  the  input,  and  the  sequence  of  reading  from  cards.  Thus,  for  example,  the 
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programmer  could  read  in  a  vector,  operate  on  it,  store  it  on  the  drum,  and 
then  atart  over  with  the  next  vector  to  be  read  from  cards,  without  concern- 
ning  himself  at  any  point  with  the  decimal  nature  of  the  numbers  on  the  cards. 

13-2  GENERAL  DECIMAL  OUTFJT  ROUTINE 

Originator  and  Coder;   Irving  H.  Eabinowitz 

This  routine  was  designed  to  allow  almost  complete  freedom  of  output 
format  without  the  necessity  of  the  programmer's  counting  or  bookkeeping. 
The  output  format  is  described  as  a  linear  function  of  two  variables,  called 
i  and  j,  and  various  numbers  specifying  number  length,  presence  or  absence 
of  signs,  etc.  The  routine  does  its  work  essentially  by  setting  up  a  corres- 
pondence between  the  given  arrangement  of  binary  data  within  the  machine  and 
the  desired  arrangement  of  data  on  cards.  To  output  a  group  of  data,  it  is 
necessary  to  specify  six  words  of  information.  These  are 

1)  The  location  and  arrangement  of  the  binary  data  on  the  drum.  This 
word  consists  of  the  coefficients  of  the  formula 

L(n. .)  =  L   +  ai  +  bj 
^  ij'    oo        " 

The  routine  can  then  use  this  word  to  find  the  binary  data  prior  to  conversion. 

2)  The  card  number  (counting  from  zero)  to  which  the  decimal  version  of 
n. .  is  to  go.  This  word  consists  of  the  coefficients  of 

C(n. .)  =  C   +  ci  +  dj 
ij'  oo        "^ 

3)  The  column  of  the  card  in  which  the  decimal  number  (including  its 
signs,  if  any,  and  its  decimal  point,  if  any)  is  to  start.  Again,  this  word 
consists  of  the  coefficients  of 

K(n. .)  =  K   +  ei  +  f j 
xj'  oo        '' 

k)     The  limits  of  i  and  j .  These  are  the  numbers  I  ,  I^ ,  J  ,  and  J^  for 

which  Iq  <  i  <  I^  and  J^  <  j  <  J^. 

5)  A  word  specifying  the  representations  of  the  plus  sign,  the  minus 
sign,  the  decimal  point;  also  the  number  of  digits  of  the  output,  and  the 
number  of  digits  in  front  of  the  decimal  point. 
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6)  Fiaally,  a  word  telling  how  many  cards  are  to  be  punched  out. 

Aa  an  example  to  illustrate  the  U3e  of  the  routine,  let  ua  consider  the 
problem  of  converting  a  5^5  matrix  from  binary  to  decimal.  Suppose  that  the 
matrix  is  stored  on  the  drum  at  locations  (08,0^4-, 00)  +  j  +51?  i.e.  the  ele- 
ments of  a  row  are  in  consecutive  locations,  and  the  various  rows  follow  each 
other  consecutively.  Since  we  would  like  to  tabulate  the  output  in  the  form 
of  a  matrix,  let  us  put  the  five  numbers  of  a  row  onto  a  single  card,  row  0 
going  to  the  first  card,  row  1  to  the  second,  etc.  Thus  we  aem   that  the  card 
number  is  simply  C(n. .)  =  i.  Since  we  have  five  numbers  to  put  on  a  card, 
we  can  break  the  card  up  into  five  fields,  each  eight  columns  wide,  starting  at 
columns  0,  8,  l6,  2k,   and  32.  These  fields  will  then  represent  a  column  j  of  the 
matrix,  and  we  can  say  that  the  starting  card  column  number  for  matrix  column  j 
is  K(n. .)  =  8j.  Since  we  have  a  5^5  matrix,  we  obviously  have  I  =  J  =0  and 

*lj'"  '  '^  oo 

I  =  J-,  =  ^.  Now  we  must  decide  on  the  format  of  the  individual  numbers.  Let 
ua  suppose  that  we  will  waat  a  blank  column  to  represent  a  plus  sign,  a  "12" 
punch  to  represent  a  minus  slgnj  that  we  do  not  wan'c  a  decimal  point,  and  that 
we  want  four  decimal  places  in  the  answer.  Then  the  fifth  word  will  contaia 
the  information  that  E(+)  =  0,   E(-)  =  "12",  E(.)  =  0,  D  ==^,  where  E(x)  stands 
for  "the  representation  of  x".  Finally,  we  will  want  to  puach  out  five  cards. 

The  six  words  will  then  be  punched  onto  cards  (in  binary)  as  follows; 

1)  00 /08,0'+, 00/00,05/00, 01        Location 

2)  00,00/00,01/00,00/00,00        Card  number 

3)  00,00/00,00/00,08/00,00        Column  number 
h)      00, 00/00, Oi)-/00, 00/00, Oi)-        Limits 

5)  00,00/00, 12/00, 00/00, Oi<-        Eepresenta+ions  and  length 

6)  00,05/00,00,00,00,00,00        Number  of  cards  to  be  punched 
The  programmer  then  inserts  this  card  into  the  code  deck  and  loads  the 

code,  assuming  that  the  binary  data  is  already  stored  on  the  drum,  and  out 
comes  the  decimal  data. 

This  is  a  particularly  simple  example  of  the  use  of  the  routine.   It  has 

provisions  in  it  for  allowing  ^he  programmer  to  write  a  "de-scaling"  routine, 

to  call  the  output  routine  from  the  drum,  to  have  the  output  routine  handle 
as  many  as  22  different  groups  of  output  cards,  and  to  have  complete 
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control  of  the  proceaa  of  output,  even  to  the  point  where  he  may,  in  a 
coded  manner,  change  the  parameters  of  the  output. 

The  code  has  a  fairly  complete  search  for  consistency  among  the 
parameters  so  that,  for  example,  the  programmer  cannot  ask  for  more  than 
forty  columns  of  information  on  a  single  card.  Stops  are  provided  at 
these  error  discoveries  with  sufficient  information  to  allow  the  pro- 
grammer to  either  find  his  error  easily  off  the  machine,  or,  if  he  recogn- 
izes it  immediately,  to  correct  it  by  simply  putting  the  correct  parameter 
word  into  the  machine  via  the  av-^cumulator. 

It  is  hoped  that  a  new  version  of  the  General  Decimal  Input  and 
General  Decimal  Output  will  be  written  in  the  near  future,  and  that  it 
will  then  no  longer  be  necessary  for  a  programmer  to  concern  himself  with 
these  processes. 
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lU.O  ASBY  -  AN  ASSEMBLY  CODE 

Originator  and  Coder;  Bryant  Tuckerman 

The  necessity  for  some  sort  of  aaseiably  code  to  be  used  ia  conjunct- 
ion with  an  automatic  computer  has  long  been  recognized,  for  such  a  code 
has  the  properties  that  not  only  is  coding  and  debugging  made  simpler,  but 
it  allows  the  building  of  a  library  of  subroutines  which  may  be  incorporated 
into  any  cede,  without  the  necessity  for  readdressing.  Furthermore,  if  it 
is  necessary  t®  change  a  code,  the  fact  that  it  is  written  in  a  language  more 
general  than  real  machine  code  ia  a  great  advantage,  since  changes  to  one 
region  of  a  code  need  not  in  general  affect  other  regions. 

The  ideal  type  of  assembly  code  is  of  course  that  one  which  takes  the 
mathematics  and  English  that  the  coder  writes  aa  his  statement  of  the  problem, 
and  turns  this  into  machine  coding  without  the  necessity  for  human  interven- 
tion. However,  such  codes  are  extremely  difficult  to  accomplish,  and,  Indeed, 
their  theory  is  not  fully  understood.  It  is  therefore  necessary  to  restrict 
ourselves  to  somewhat  less  sophisticated  types  of  assembly  routines.  Among 
other  considerations  that  must  be  applied  are  those  of  time  and  space.  Since 
machine  time  is  a  valuable  commodity,  the  assembling  of  a  code  written  in 
some  external  language  should  be  only  a  very  small  proportion  of  the  running 
time  of  the  assembled  code.  Furthermore,  siace  space  is  at  a  premium,  the 
secondary  storage  of  the  machine  being  finite  (and  for  eoms  applications, 
quite  small),  it  is  desirable  that  the  input,  or  unassembled  version  of 
the  code  not  be  much  larger  than  its  final  version.  In  order  to  reach  these 
objectives,  certain  compromises  between  efficiency  of  machine  use  and  ease 
of  coding  must  be  made. 

ik . 1  Format  of  Subroutines 

ASBY  assumes  that  all  coding  consists  of  subroutines  written  in  a  special 
language  called  Format  C  which  is  almost  machine  language.   In  general,  address- 
es are  relative,  i.e.  they  consist  of  a  tag  identifying  the  subroutine  and  a 
number  which  indicates  the  address  within  that  subroutine.  Addresses  refer- 
ring to  words  within  the  present  subroutine  are  all  larger  than  (l6,00). 
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i.e.  all  subroutines  are  written  as  if  they  were  to  be  operated  from  a 
block  of  consecutive  memory  locations  starting  from  (l6,00).  Furthermore 
no  subroutine  may  be  longer  than  8  lines  (256  words).  Addresses  refer- 
ring to  other  subroutines  are  smaller  than  (l6,00)  and  are  subdivided  into 
Ik   categories,  according  to  the  first  pentad  of  the  address  with  0,1  and 
2  being  treated  as  a  single  category.  The  addresses  beginning  with  0,1,2 
are  treated  as  absolute  addresses  and  are  never  modified  by  ASBY-  The 
other  thirteen  categories  of  addresses  begin  with  the  integers  3  through 
15,  which  may  be  considered  13  tags  identifying  13  other  subroutines.  The 
rest  of  the  address  is  a  pentad  (O  through  31) .  thereby  limiting  all  refer- 
ences to  other  subroutines  to  the  first  32  instructions  in  these  subroutines. 
This  is  no  real  restriction  for  logical  operation  since  it  is  difficult  to 
imagine  a  subroutine  which  will  require  more  than  32  different  types  of  ex- 
ternal references  to  it;  usually  2  or  3  will  suffice.  Data,  on  the  other 
hand,  must  also  be  organized  into  our  subroutine  format  if  they  are  to  be 
assembled  together  with  the  code  and  here  the  limitation  of  only  32  distinct 
references  may  require  some  indirect  addressing  system  for  blocks  of  data. 
Data  subroutines  can,  of  course,  be  256  words  long. 

While  any  one  subroutine  will  probably  not  need  to  refer  to  more  than 
13  other  subroutines,  nevertherless  we  shall  very  quickly  acquire  a  library 
of  considerably  more  than  13  subroutines  with  many  mutual  cross  references. 
Thus  it  is  not  practical  to  limit  a  tag  (say  tag  5)  to  refer  always  to  the 
same  subroutine  (say  SINE  of  X)  whenever  it  is  used.  Format  C  therefore 
provides  a  tag-defining  group  of  words  which  sits  at  the  front  end  of  the 
subroutine  so  long  as  the  subroutine  remains  in  the  library,  but  which  is 
removed  from  the  subroutine  after  serving  its  purpose  during  assembly.  This 
tag-defining  group  consists  of; 

1)  a  single  parameter -word,  PARW,  having  ones  in  those  binary 
stages  corresponding  to  the  tags  used  in  this  subroutine,  and  zeros 
elsewhere.   (if  we  use  tags  3,^  and  7  then  PABW  is  OOOllOOlO. • -O.) 

2)  a  sequence  of  IDENtifing  words,  one  for  each  tag  used, 
arranged  in  increasing  numerical  order QOfiBl&et*9g8.  (In  our  example 
above,  PAEW  would  be  followed  immediately  by  three  IDEN  words,  the 


first  bearing  the  identificatio*  of  that  subroutine  to  which  we  are 
referring  hy  the  tag  3,  the  next  Identifying  the  subroutine  we  are 
invoking  with  tag  h,  and  the  third  identifying  the  subroutine  we 
refer  to  by  the  tag  7«) 

In  addition  to  these  tag-defining  words,  a  one-word  description, 
DSCE,   of  the  current  subroutine  also  appears   in  the  front-end  of  the 
library  form  of  all  Foriaat  C  subroutines.     This  word  specifies  the  Williams 
assembled  length,   dw,   of  the  subroutine;  the  library  (Format  C)    length,   d£, 
of  the  subroutine;  and  the  number  of  constants,   c,  at  the  end  of  the  sub- 
routine which  are  to  be  left  unmodified.     Each  of  these  discriptive  numbers 
occupies  10  bits.     We  may  sumEiarize  this  form  by 

DSCE  =  (dw;   di  ;   c;  0) 

All  this  front-end  iaformatioa  of  ASBY  is  placed  between  words  (l6,00) 
and  (l6,01)  and  none  of  it  acquires  any  location  address  since  it  will  ul- 
timately be  removed.  Word  (1^^,00)  is  the  full  word  IDENtification  of  the 
subroutine  and  the  actual  code  begins  in  (l6,01).  Thus  we  have  our  complete 
Format  C  outline; 

(16.00)  IDEN  (  ) 
DSCR  (  dw;  dl  ;  c;  0  ) 
PAEW  (  ) 
IDEN^  (  ) 

(16.01)  code   (  ) 

(     )  B§^e   (  ) 

-SUM   (  ) 

The  word  designated  by-SUM  is  for  checking  purposes,  and  contains  the 
negative  of  the  sum  (modulo  2)  of  all  the  other  words  of  the  subroutine.  No 
gaps  in  the  code  within  a  subroutine  are  permitted  and  all  subroutines  are 
punched  in  straight  binary,  12  words  per  card  except  possibly  at  each  end 
where  aa  entire  card  may  not  be  needed.  Routines  may,  of  course,  be  punched 
as  decimal  pentads  and  then  converted  and  punched  out  in  binary  by  using 


existing  codes.  The  total  number  of  noa-Willlaa  nemory  words  is  3  +  t 
vkere  t  is  the  nuaber  of  tags  used,  hence 

d£  =  dw  +  3  +  t 

A  special  feature  of  ASBY  allows  us  to  define  fields  of  variables  of 
arbitrary  but  given  length,  dw,  by  subroutines  occupying  only  d£  =  k   library 
words ,  in  the  form 

(16,00)   IDEN  (  ) 

DSCE  (  dw;  Uj  0;0   ) 

PAEW  (      0  ) 

-STJM  (  ) 

This  is  assembled  as  dw  words,  the  first  being  IDEN,  the  rest  0  or  trash 
according  as  the  drum- image  of  Williams  was  initially  cleared  or  not. 

(More  generally,  any  number  of  irrelevant  or  zero  words  in  the  end  of 
a  subroutine  may  be  omitted  from  the  library  form  to  become  0  or  trash  upon 
assembly,  by  proper  choice  of  di'  <  dw  +  3  +  t.) 

lU.2  Directory  Format  for  ASBY 

Im  addition  to  his  library  of  subroutines  the  customer  must  also  supply 
a  Directory  which  must  include  at  least  the  IDEN  of  the  first  routine  to  be 
obeyedo   In  general  no  furthftr  information  need  be  supplied,  'since  ASBY  will 
look  at  thia  routine,  will  see  what  others  are  called  for,  will  examine  them 
in  turn,  etc.  until  it  has  constructed  the  complete  directory  during  assembly. 


*)  A  non-zero  starting  Williams  address  should,  however,  be  supplied  by  the 
customer  for  the  location  of  the  first  subroutine  as  w  of  its  INFO  word. 
Otherwise  the  first  routine  will  go  into  (00,00),  thereby  violating  the 
convention  that  we  will  keep  the  bottom  12  words  of  the  Williams  memory 
free  for  small  service  routines. 
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Tlie  customer  may,  kcwover,  override  part  of  this  process  by  supplyiag 
tke  IDEN  aad  assembled  locatiom  of  any  of  the  other  subroutimes  he  wishes 
to  neatio*.   (This  facility  may  turm  out  to  be  useful  duriag  debugging, 
simce  he  will  thea  kaow  that  certain  crucial  subroutines  begia  ia  some 
easily  remembered  locatioas.   It  may  also  be  used  to  cause  the  loadiag  of  auxil- 
iary subroutiaes  aot  directly  aeeded  by  the  code  but  which  would  be  valuable 
to  have  available  ia  the  Williams  memory  duriag  debuggiag.) 

The  final  form  of  the  Directory  coasists  of  pairs  of  words 


IDEN^ 

( 

(    WJ    d£    J 

U5     / 

) 
) 

HJENg 
iHJ-'Og 

( 

(  wj   d^  1 

UJ      / 

) 

) 

where  w  (10  bits^  =  Williams  address  of  the  asseiibled  subroutiae 
dl   (10  bits^  =  library  leagth  of  the  subroutine 
u  (5  bits)  =  unimportant  and  irrelevant 
^  (15  bits)=:  drum  address  of  the  subroutine 

If  the  customer  wishes  to  fix  the  location  of  a  subroutiae,  he  need  only 
list  its  IDEM  and  then  supply  its  w  (the  other  information  will  be  constructed 
correctly  by  ASBY  even  if  supplied  here  incorrectly.) 

In  the  absence  of  w,  ASBY  will  supply  it  as  the  next  conaecutuve  location. 
However,  at  the  conclusion  of  assembly  the  completed  Directory  is  punched  out 
and  control  passes  to  the  RELOAD  code. 


II.      SOLUTION     OF       PEOBLEMS 


20.     General  Remarks 

During  the  period  covered  by  thie  report  the  emphasis  has  not 
primarily     been  placed  on  "production"  but  ratter  on  exploring  the 
potentialities  of  the  computer  as  a  research  tocl  in  a  variety  of 
problems  and  applications,  while  our  engineering  group  was  engaged 
in  improving,  hence  changing  our  equ.ipmento     As  a  result,  no  syste- 
matic effort  was  made  at  the  tiiae  tc  bxiild  up  a  library  of  subroutines, 
or  to  standardise  coding  procsdures,  and  the  number  of  problems  solved 
is   comparatively  small. 

It  would  not  be  easy  to  arrange  the  probleias  treated  on  our  machine 
into  clear-cut  groups i  three  such  groups  can  be  distinguished,  however, 
and  we  shall  deal  with  them  separately. 

The  first  should  be  headed  ME1S0E0L0GY.     Dr.   Charney's  group  used 
about  one  third  of  the  total  machine  tin»  available  for  operation  during 
the  period  covered  by  this  report.     Its  problems  and  results  will  be  re- 
ported under  the  terms  of  Contract  Nonr  1358-(02)    (between  the  Institute 
for  Advanced  Study  and  the  Office  of  Naval  Research) ,  by  which  the  meteo- 
rological research  and  machine  coding  were  supported. 

The  second  goup  of  problems  deals  with  ASTROPHYSICS.     All  of  this 
work  has  been  directed  by  Prof.  Martin  Schwarzschild   (Princeton  University) , 
either  directly  or  via  his  co-workers.     These  problems  are  discussed  in 
chapter  21,  to  which  Prof.   Schwarzschild  was  kind  enough  to  write  an  in- 
troduction. 

The  problems   in  the  third  group  are  only  loosely  related  to  each 
other  by  virtue  of  their  common  subject  matter,  viz.  ATOMIC  AND  NUCLEAR 


PHYSICS.  No  two  problems  have  originated  with  the  same  individual.  It 
seemed  advisable,  however,  to  \inite  them  in  one  chapter  (22)  so  that  the 
reader  interested  in  physics  can  find  them  more  easily. 

The  last  chapter  (23)  is  nothing  but  a  collection  of  those  problems 
which  do  not  fit  into  any  of  the  other  groups.  Size  and  importance  of  these 
problems  are  as  varied  as  their  subject  matter.  The  Historical  Ephemeria 
has  taken  close  to  one  man-year  of  coding  time  (nearing  completion  at  the 
close  of  this  report)  and  will  be  of  lasting  value.  The  remaining  problems 
are  considerably  shorter,  and  a  small  number  of  minor  codes,  without  any 
special  features,  have  been  omitted.  They  have  been  helpful,  however,  in 
demonstrating  the  importance  of  library  subroutines  and  determining  their 
specifications . 
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21.  ASTROPHYSICS 


Astrophysics  differs  from  other  physical  sciences  by  the  circumstance 
that  its  objects  of  study  can  be  observed  but  cannot  be  experimented  with. 
This  impossibility  of  physical  experimentation  can  be  compensated  for  to  a 
remarkable  degree  by  numerical  experimentation.   It  is  for  this  reason  that 
numerical  research  plays  such  an  unusually  important  role  in  astrophysics. 


The  problems  here  investigated  are  of  two  types;  the  first  type  refers 

th< 
2) 


to  the  undisturbed  internal  structure  and  evolution  of  stars,  '  whereas  the 


second  type  refers  to  perturbations  of  stars  such  as  stellar  pulsations 
or  hydromagnetic  waves. '^  The  first  type  leads  to  highly  non-linear,  high 
order  eigenvalue  problems.   In  them  the  main  difficulty  consists  in  finding 
efficient  methods  of  determining  the  eigenvalues  by  trial  and  error.  The 
second  type  leads  to  partial  differential  equations  with  as  much  as  three 
independent  dimensions  (two  in  space  and  one  in  time).  Here  the  main  dif- 
ficulty consists  of  finding  efficient  methods  for  the  numerical  intergra- 
tion  that  are  free  from  numerical  instabilities. 


1)  cf.  §§  (21.10),  (21.43) 

2)  cf,  §§   (21.20),  (21.50) 

3)  cf.  f   (21.30) 
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21.10  NUMERICAL  EXPERIMEWTATION  ON  STELLAE  EVOLUTION 

Originator  and  Analsrat:  M.  Schwarza child 
Coder;  Mrs.  H.  Selberg 

Recent  progress  in  nuclear  physics  has  made  it  possible  to  determiae 
fairly  accurately  the  reaction  rates,  as  functioms  of  temperature  and 
density,  of  those  nuclear  processes  which  provide  the  main  energy  sources 
in  the  stellar  interior.  These  reaction  rates  have  provided  the  last  link 
necessary  to  formulate  the  problem  of  stellar  structure  and  evolution  in 
a  unique  manner.  Consequently,  it  is  now  possible,  in  principle,  to  derive 
all  physical  characteristics,  such  as,  for  example,  central  temperature  and 
total  luminosity,  for  a  star  of  given  mass  and  initial  composition  as  a  func-  u 
tion  of  time  throughout  the  star's  evolution. 

In  practice,  however,  the  solution  of  this  problem  is  complicated, 
both  as  regards  the  methods  of  numerical  analysis  to  be  applied,  and  as 
regards  the  effective  carrying  through  of  the  large  amount  of  numerical  cal- 
culations. The  difficulties  arise  largely  from  the  fact  that  in  the  space 
co-ordinate  (from  the  center  to  the  surface  of  the  star)  the  physical  condi- 
tions present  a  highly  non-linear  system  of  differential  equations  of  the 
fourth  order  with  two  boundary  conditions  at  the  center  and  two  at  the  surface, 
but  also  from  the  fact  that  the  energy  release  by  gravitational  contraction 
gives  the  problem  in  the  time  co-ordiaate  the  character  of  a  heat  coaduction 
problem  with  its  well-known  danger  of  numerical  instability. 

The  numerical  solution  of  problems  in  the  theory  of  stellar  evolution 
is  clearly  of  importance  to  theoretical  astrophysics.  Simultaneously,  the 
development  of  methods  for  the  solution  of  this  problem  appears  likely  to  be 
valuable  for  the  solution  of  a  large  class  of  non-linear  high  order  boundary 
value  problems. 

21.11  Qualitative  Description  of  Stellar  Evolution 

Consider  a  star  of  a  given  mass  and  a  given  initial  chemical  composition. 
Assume  that  the  matter  of  which  the  star  consists  was  originally  well  mixed  so 
that  the  star  starts  out  chemically  homogeneoiis.  After  a  pre-stellar  contraction 
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phase  a  star  will  settle  into  its  inital  equilibrium  configuratioa  in  which 
all  the  hydrostatic  and  thermal  equilibrium  conditiona  are  fulfilled.   In 
particular  the  internal  temperatures  will  be  just  right  so  that  the  hydrogen 
burning  produces  energy  at  a  rate  exactly  compensating  the  losses  by  surface 
radiation. 

The  star  would  remain  indefinitely  in  this  particular  equilibriua 
configuration  if  it  were  not  for  the  circumstance  that  in  time  the  nuclear 
burning  alters  the  chemical  composition  of  the  deep  interior,  which  in  turn 
causes  changes  in  the  equilibrium  configuration.  The  investigation  of  stellar 
evolution  aims  to  derive  the  sequence  of  configurations  through  which  a  star 
evolves  and,  in  particular,  to  compute  the  changes  of  the  luminosity  and  the 
radius  of  a  star  during  its  evolution.  A  comparison  of  these  two  quantities 
with  observation  provides  a  check  on  theoretical  developments. 

A  previous  investigation  ^^has  suggested  the  following  sequence  of 
evolutionary  events  for  a  star  with  a  mass  similar  to  that  of  the  sun.  In 
the  initial  phases  more  and  more  hydrogen  in  the  core  is  transmuted  into 
helium.  After  a  while  the  hydrogen  in  the  core  is  exhausted  and  the  hydrogen 
burning  continues  in  a  shell  further  out.  The  helium  core  becomes  isothermal. 
It  steadily  grows  in  mass  but  shrinks  in  size  while  simultaneously  the  envelope 
of  the  star  becomes  more  and  more  expanded.  The  contraction  of  the  core 
produces  steadily  increasing  degeneracy  of  the  gases  near  the  center  while 
the  growth  of  the  envelope  causes  a  convective  zone  below  the  stellar  surface 
which  steadily  increases  in  depth. 

During  these  evolutionary  phases,  both  the  luminosity  and  the  radius 
of  the  star  continuously  increase,  i.e.  the  star  becomes  a  red  giant. 

As  the  luminosity  increases,  the  rate  of  nuclear  burning  as  well  as 
the  rate  of  the  core  contraction  will  increase.  Eventually  the  contraction 
rate  will  be  so  high  that  the  normal  thermal  equilibrium  condition  does  not 
apply  any  fxirther.  The  contraction  will  heat  the  center  of  the  core  by 


*)  F.  Hoyle  and  M.  Schwarzschild,  Ap.  J.,  Suppl.  No.  13,(1955) 
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compressio*  -  quite  Irrespective  of  the  lack  of  kydroge*  fuel  for  nuclear 
buraimg  -  until  it  reackes  a  temperature  sufficient  for  helium  burning  to 
commeace.  This  appears  to  be  the  critical  phase  at  which  the  star  has  reach- 
ed the  top  of  the  observed  red  giaat  sequence.  Whea  the  helixim  buraiag  starts, 
it  will  contribute  to  the  rise  in  the  central  temperature.  This  rise  will  aot, 
as  la  the  aoa-degeaerate  case,  be  checked  by  an  immediate  expaasioa  of  the  core 
siace  the  now  dominant  degenerate  pressure  is  independent  of  the  temperature. 
Thus  an  unstable  situation  may  arise  for  a  short  evolutionary  phase.  During 
this  phase,  it  appears  possible  that  the  temperature  gradients  may  become  too 
high  to  be  staolr  against  convsctioa  so  that  coavectioa  may  set  in  throughout 
most  of  the  star,  which  will  mix  the  helium  of  the  core  with  the  hydrogen  of 
the  envelope.  Such  a  thorough  mixing  would  obviously  have  important  consequences 
for  the  interpretation  of  the  subsequent  evolutionary  phases. 

So  far,  the  acc^orate  sequence  of  events  d-uring  the  critical  phases  at 
the  top  of  the  red  giant  sequence  have  act  been  determined  siace  the  onset  of 
helium  burning  and  the  active  influence  of  contractive  heating  make  the  computa- 
tions difficult.   It  is  the  purpose  of  this  project  to  investigate  these  fast 
and  critical  evolutionary  phases. 

21.12  The  Basic  Equations  and  Boundary  Conditions 

The  basic  conditions  throughout  the  stellar  interior  are  governed  by 
three  differential  equations.  The  first  one  expresses  the  fact  that  hydro- 
static equilibrium  holds,  i.e.  that  the  force  produced  by  the  gas  pressure 
exactly  balances  gravity.  The  second  one  expresses  the  relation  between  the 
energy  flux  and  the  temperature  gradient.  Here  three  cases  have  to  be  distin- 
guished, depending  upon  whether  the  flux  is  transported  by  radiation,  by 
convection  or  by  electron  conduction,  the  latter  being  the  case  when  the  gas 
is  degenerate.  The  third  differential  equation  expresses  the  conservation  of 
energy;  the  divergence  of  the  energy  flux  must  be  equal  to  the  energy  genera- 
tion by  nuclear  processes  to  which,  in  the  phases  of  fast  contraction,  the 
energy  freed  by  non-adiabatic  contraction  is  to  be  added. 
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This  system  of  differential  equations  contaias  three  auxiliary 
variables s  the  density,  the  absorption  coefficient,  and  the  rate  of 
energy  geaeration,  which  are  related  to  the  main  variables  by  rather 
complicated  equations. 

The  density  is  given  as  a  function  of  temperature  and  density  and 
of  the  abundances  of  helium  and  hydrogen.  There  are  two  such  equations  of 
state,  one  for  the  non-degenerate  and  one  for  the  degenerate  case,  and  a 
rather  involved  criterion  must  be  used  to  decide  which  equation  to  apply. 
The  switch  ia  presently  made  abruptly,  while  accurately  an  intermediate 
equation  for  incipient  degeneracy  should  be  used  --  an  accuracy  which  is 
hardly  required  here. 

The  absorption  coefficient  is  given  as  a  sum  of  two  terms.  The  first 
represents  electron  scattering  and  does  not  depend  on  the  temperature;  the 
second  terra  conbines  free -free  transitions  for  hydrogen  and  helium  with 
bound -free  transitions  for  the  heavier  elements  and  is  proportional  to  t"   . 
The  simple  addition  of  the  terras  is  an  approximation,  again  of  sufficient 
accuracy  for  the  present  purpose. 

The  rate  of  energy  generation  is  again  given  as  a  sum  of  two  terms 
representing  the  energy  liberated  by  burning  hydrogen  and  helium,  respectively. 
Probably  these  two  terms  will  never  contribute  simultaneously  in  the  same 
volume  since  presumably  the  hydrogen  will  be  completely  exhs-uated  before  the 
temperature  reaches  the  threshold  for  helium  burning  to  set  in.  Both  terms 
depend,  of  course,  not  only  upon  the  temperature  but  also  on  the  respective 
abundances  of  hydrogen  and  helium  which,  however,  must  be  replaced  by  mean 
values  for  the  convective  region  where  perfect  mixing  may  be  assumed. 

21.13  Boundary  Conditions  and  Automatic  Adjustment  of  Boundary  Values 

To  complete  the  definition  of  the  problem  the  boundary  conditions 
at  the  center  and  at  the  surface  of  the  star  have  to  be  assembled. 

Our  four  dependent  variables  are;  the  pressure  P(r) ,  the  temperature 
T(r),  the  mass  M  contained  in  the  sphere  of  radius  r,  and  the  luminosity 
L  ,  i.e.  the  amount  of  energy  radiated  per  unit  time  through  the  surface 
of  that  sphere.  The  boundary  conditions  at  the  center,  therefore,  are 
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M  =  L  =  0  while  P  amd  T  are  aot  kaovfa  a  priori.  Tke  other  two  bouadary 
comditioHS  follow  from  the  equilibrium  conditions  at  the  Surface,  i.e.  in 
the  photosphere  of  the  star;  the  first  states  the  relation  between  the 
surface  temperature  and  the  total  luminosity  and  the  second  insures  that 
the  pressure  in  the  photosphere  provides  the  appropriate  optical  depth 
of  this  layer. 

To  introduce  these  two  surface  conditions  in  their  accurate  form  would 
involve  appreciable  practical  difficulties  because  it  would  necessitate, 
first,  the  addition  of  the  whole  set  of  photoapheric  equations  including 
partial  ionization  of  hydrogen,  and  second,  the  handling  of  the  large  ranges 
in  temperature  and  pressure  from  the  center  of  the  star  to  its  photosphere. 
These  complications  are  here  avoided  by  the  following  approximating  procedure* ■ 
The  outermost  layers  of  the  star  are  cut  off  at  the  point  where  the  temperature 
reaches  exactly  one  million  degrees.  This  cutoff  reduces  the  total  mass  only 
by  an  entirely  negligible  amount.  It  reduces  the  radius  by  a  somewhat  larger 
percentage,  which,  however,  la  still  not  serious  and  can  be  corrected  for, 
if  need  be,  at  the  end  of  the  calculations.  The  most  serious  consequence 
of  this  procedure  is  that  the  pressure  at  the  cutoff  point  has  to  be  estimated 
from  previous  tentative  calculations. 

If  all  of  the  four  boundary  conditions  were  given  at  the  same  end  of 
our  integration  interval,  a  single  numerical  integration  would  yield  the  con- 
figuration at  each  definite  phase «  But,  unfortunately,  this  is  not  the  caaej 
we  have  just  seen  that  we  have  two  conditions  at  each  end.  Therefore,  a 
trial  and  error  method  must  be  \iaed  to  determine  the  initial  values.  For 
example,  we  could  assume  certain  estimated  values  for  P  and  T  at  the  center 
and  integrate  out  to  the  surface.  We  could  then  vary  P(o)  and  T(o)  in  some 
systematic  way  until  the  computed  values  would  fulfill  the  boundary  conditions 
at  the  surface.  Or  we  could  start  the  integrations  with  estimated  values  at 
the  surface  (fulfilling,  of  course,  the  two  boundary  conditions)  and  Integrate 
in  until,  by  a  proper  choice  of  the  two  remaining  free  parameters  at  the  surface, 

M  and  L  would  come  out  close  enough  to  zero.   --It  turned  out,  however,  that 

r      r 
neither  of  these  schemes  is  applicable  due  to  the  inherent  instability  of  our 

system  of  differential  equations  as  soon  as  the  values  of  the  dependent  variables 


21.11^.0 

depart  from  tlie  pkysically  correct  combinatioa.  However,  a  combiaatio* 
of  both  methods  was  successful.   It  basically  consists  of  iategrating 
from  both  eads  toweurda  some  iatermediate  "fitting  poiat"  for  which  the 
variables  caa  be  estimated  fairly  well.   The  free  parameters  at  either 
ead  are  first  varied  uatil  the  values  at  the  fitting  poiat  are  aear 
the  estimated  values.   The  final  fittiag  is  thea  made  by  liaear  iater- 
polatioa. 

It  is  aot  possible  to  describe  the  details  of  this  fittiag  pro- 
cedure here,  but  we  should  like  to  meation  that  Prof.  Schwarzschild  and 
Mrs.  Selberg  succeeded  ia  plaaniag  aad  codiag  this  procedure  so  that  ia 
spite  of  the  highly  aoa-liaear  character  of  this  "eigeavalue"  problem 
It  ia  carried  out  by  the  Computer  without  human  iaterference  except  for 
feeding  the  initial  guesses  for  the  variables  at  the  center,  the  fitting 
point  and  the  surface  (i.e.  the  one  million  degree  point). 

21.  l**-  General  Computational  Procedure 

The  numerical  integration  of  our  system  of  differential  equations  is 
hampered  by  its  high  degree  of  non-liaearity  and  the  very  wide  range  of  al- 
most all  variables.  Furthermore,  the  whole  system  including  the  various 
supplementary  equations  for  computing  the  opacity,  energy  generation,  etc. 
and  the  distinction  between  the  radiative,  convective  aad  degeaerate  cases 
is  certaialy  the  most  complicated  sjrstem  ever  treated  oa  our  computer.   It 
is  impossible,  therefore,  to  give  a  full  account  of  the  computatioaal  pro- 
cedure ia  this  short  report.  However,  a  few  remarks  may  be  appropriate  here. 

To  start  with,  the  spatial  iadepeadent  variable  is  changed  from  r 

to  M  .  This  is  essentially  a  transformation  from  Eulerian  to  Lagrangian 

r 
coordinates,  which  makes  it  much  easier  to  follow  the  changes  in  chemical 

composition  by  nuclear  transmutations  in  each  mass  element. 

Next,  let  us  consider  the  derivatives  with  respect  to  time.  Such 
derivatives  occur  in  the  law  of  energy  conservation  and  in  the  law  for  trans- 
mutation rates.  The  energy  conservation  equation  has  the  character  of  a 
heat  conduction  equation.   Therefore,  its  explicit  use  for  extrapolation 
in  tine  in  a  step-wise  numerical  integration  leads  to  a  well  known  instability 
except  if  very  small  time  steps  are  used,  which  would  be  impracticable  in 
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tke  preseat  applicatiom.   Co»sequemtly,  tkis  equatioia  must  be  used  im  am 
"implicit"  naaiier,  i.e.  wke»  aa  evolutioaary  time  step  leadiag  from  pre- 
vious stellar  model  to  tke  next  stellar  model  is  to  be  bridged,  tke  equatioa 
of  emergy  coaservatioa  kas  to  be  applied  ia  a  symmetrical  maaaer  botk  at  tke 
previous  time  aad  tke  aew  time.   Oa  tke  otker  kaad,  tke  derivatives  witk 
respect  to  time  occurriag  ia  tke  equatioas  goveraiag  tke  traasmutatioa  rates 
may  be  used  ia  aa  explicit  maaaer  for  extrapolatioa  ia  time,  ao  reasoa 
beiag  kaowa  to  siispect  kere  a  aumerical  iastability.  Tke  most  simple  liaear 
extrapolatioa  formula  would  be  aufficieat  for  tke  accuracy  required  kere, 
but  tt  may  lead,  oa  occasioa,  to  erroaeous  aegative  values  for  tke  kydrogea 
or  keliura  abuadaace.  A  suitable  correctioa  was  used  to  avoid  tkis  coaditioa. 

Altogetker  tke  kaadliag  of  oae  evolutioaary  time  step  is  tkus  reduced 
to  tke  followiag  computatioa;  after  tke  derivatioa  of  tke  coafiguratioa  for 
oae  defiaite  time  kas  beea  completed,  tke  coatractioa  aad  tke  aew  kydrogea 
aad  kelium  abuadaacer are  computed  for  every  poiat  ia  tke  star.  Witk  tkeae 
tkree  fuactioas  tke  time  laterval  is  bridged  aad  tke  coafiguratioa  at  tke 
aext  ^kaae  caa  be  derived. 

To  adjust  all  tke  equatioas  to  tke  requiremeats  of  tke  electroalc 
computer,  tke  pkysical  variables  kad  to  be  scaled  appropriately.  Most  of 
tkem  were  scaled  liaear ly  i.e.  multiplied  by  aa  appropriate  power  of  2 
or /aad  10^  ia  additioa  tke  Aenaitj  §   aad  tke  pressure  P  were  replaced  by 
tkeir  tkird  aad  fourtk  roots,  respectively,  ia  order  to  reduce  tke  raage 
of  tke  variables  kaadled  by  tke  computer.  Furtkermore,  tke  distaace  r 
from  tke  ceater  was  replaced  by  tke  variable  k  =  l/r,  wkick  kas  its  kigkest 
aumber  of  sigaiAcaat  figures  ia  tke  importaat  ceatral  regioa  of  tke  star 
ratker  tkaa  ia  tke  less  importaat  eavelope.   Tke  time  step  was  scaled  ia 
suck  a  way  as  to  iasure  approximate  coastaacy  of  tke  amouat  of  kydrogea  or 
keliua  buraed  per  time  step. 

Tke  iatroductioa  of  tke  variable  k  =  l/r  leads  to  a  (aoaHBsseatial) 
siagularity  at  tke  center  of  tke  star.  Tkerefore,  a  simple  aaalytlc  approxima- 
tioa  is  used  to  start  tke  aumerical  iategratioa  at  some  small,  but  fiaite 
radius.  From  tkereoa  tke  leagtk  of  eack  step  is  computed  by  tke  mackiae  from 
tke  requiremeat  tkat  ao  variable  must  ckaage  by  more  tkaa  some  fixed  fractioa 
of  its  value. 
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21.15  Machine  Utilization 

The  length  and  complexity  of  the  code  for  the  present  computation 
confronted  ua  with  several  problems  which  are  not  normally  encountered  in 
smaller  calculations.  We  shall  discuss  some  of  them  in  this  paragraph. 

One  of  the  most  difficult  ones  was  the  problem  of  storage.  Since 
it  was  not  possible  to  accommodate  the  entire  code,  including  all  the 
data  to  be  stored  for  each  model,  in  the  Williams  memory  {192h   words), 
several  subroutines  were  stored  on  the  magnetic  drum  and  called  into  the 
Williams  memory  whenever  needed.  This  procedure  proved  to  be  not  only 
very  time  consuming  --  especially  with  the  old  drum  which  had  an  access 
time  of  nearly  200  msec.  --  but  it  also  would  have  required  a  reliability 
of  the  magnetic  drum  and  its  amplifying  circuitry  that  could  not  be  achiev- 
ed with  the  old  driom.  Both  conaideratioixs  hold,  to  a  lesser  degree^  even 
for  the  new  drum.   It  was  tried,  therefore,  to  reduce  the  number  of  drum 
references  by  accommodating  in  the  Williams  memory  at  least  those  sub- 
routines which  are  most  frequently  needed,  i.eo  the  ones  accurring  in 
every  integration  step. 

When  the  new  drum  became  available,  a  restart  procedure  was  added 
to  the  old  code.  Storing  the  code  and  the  data  (about  3000  words)  at  the 
end  of  each  full  integration  (i.e.  from  center  to  surface),  it  was  possible, 
in  case  of  obvious  machine  trouble,  to  resume  the  computation  at  the  end  of 
the  previous  integrationo  In  addition,  an  optional  automatic  duplication 
mechanism  was  recently  Incorporated  in  our  code. 

As  mentioned  in  the  previous  paragraph,  scaling  presented  a  considerable 
problem.  In  fact,  a  fixed  point  machine  is  not  very  well  suited  to  this  type 
of  problem,  in  which  the  variables  vary  over  a  very  wide  range.   It  requires 
careful  and  accurate  scaling  in  order  to  preserve  a  sufficient  degree  of 
accuracy;  moreover,  since  the  ranges  of  the  variables  may  differ  from  model 
to  model,  the  scaling  has  to  be  changed  accordingly.  Two  pogsibilities  were 
investigated  to  solve  this  difficulty;   (i)  the  use  of  coded  floating  point 
techniques,  and  (11)  the  use  of  logarithmic  variables.  Of  the  two,  the  latter 
seems  more  adequate  since  our  formulae  contain  a  considerable  number  of  Integer 
and  fr-actlonal  powers.  It  is  planned,  therefore, 
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to  rewrite  the  code  If  the  subroutines  for  log  x,  exp  t: ,   and  log  (x  +  y) 
can  he  made  fast  enough  and  the  Williams  memory  can  accommodate  them 
during  a  vhole  integration  step.  This  would  also  solve  the  problem  of 
read-around  in  the  Williams  memory,  which  arose  in  the  evaluation  of 
fractional  powers  and  could  only  be  overcome  by  an  artificial  elow-dowa. 

So  far  we  have  computed  nine  models.  The  average  computing  time  for 
one  model  is  3+1/2  hours.  Actually,  many  more  machine  hours  were  spent 
on  debugging  and  on  experimentation  with  new  computational  methods,  and 
a  good  many  were  lost  due  to  machine  errors. 

21.16  Besults  and  Outlook 

With  the  procedure  described  above  a  sequence  of  nine  models  has  so 
far  been  computed  which,  for  a  specific  star,  covers  a  range  of  advanced 
evolutionary  phases  in  which  the  star  has  red  giant  character,  in  satis- 
factory agreement  with  astronomical  observations.  The  results  are  in 

1) 
general  agreement  with  previous  investigations,  '  but  they  indicate  that 

the  contraction  term  in  the  third  basic  equation  plays  a  more  important 

role  for  the  structure  of  the  advanced  evolutionary  phases  than  the  earlier 

rough  estimates  had  indicated. 

It  is  planned  to  continue  the  present  evolutionary  model  sequence 
up  to  and  beyond  the  phase  when  helium  burning  sets  in,  and  to  carry  through 
other  evolutionary  sequences  for  stars  of  different  masses  and  initial 
compositions,  for  the  purpose  of  obtaining  a  secure  basis  fcJr  the  theory 
of  those  most  advanced  evolutionary  phases  in  which  the  formation  of  the 
heavy  elements  is  expected  to  occur. 


1)  cf.  Numerical  Integrations  for  the  Stellar  Interior,  Hirm,  R.  and 

Schwarzschild,  M-^  Astrophysical  J.,  121  (1955)  and  Astrophysical 
J.,  Suppl.  Series  Vol.  I,  p.  319-^30"Tl955) ,  where  numerous  further 
references  are  given. 
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It  is  alao  planned  to  prepare  technical  reports  and/or  papers 
for  afcientific  journal  'whenever  an  important  phase  of  this  problem 
has  been  finished. 

The  supporting  contracts  for  this  vork  are:  Contract  No*  Da-36-03'*-- 
OBIi-l6h6   for  machine  operation.  Contract  Uc  Nonr  1358- (03)  for  all  coding 
and,  while  this  report  is  being  written.  Contract  JJonr  1358- (O^^)  for  further 
use  of  the  machine.  (January  1,  1957  -  June  30,  1956.) 


1)  presumably  the  Astrophysical  Journal  and/or  its  Supplementary  Series. 
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21.20  ADIABATIC  PULSATIONS  OF  AN  QRIGIUALLY  ISOTHERMAL  ATMOSPHERE 

Originator  and  Coder:  D.  A.  Lautmann  ' 

The  periodic  variation  of  brightness  of  certain  stars,  named  Cepheids 
after  the  most  famous  of  them,  viz.   Cephei,  is  one  of  the  most  intriguing 

problems  of  astrophysics.  The  so-called  pulsation  theory  attempts  to  explain 

2) 
these  variations  on  the  basis  of  hydrodynamic  pulsations.  '  The  theory  of  these 

pulsations  is  fairly  well  understood  for  the  linear  case,  i.e.  as  long  as  the 
amplitudes  of  the  oscills-tions  are  so  small  that  quadratic  terms  can  be  neglect- 
ed. In  the  deep  interior  of  a  star,  where  energy  considerations  demand  small 
amplitudes,  the  highly  developed  linear  theory  of  adiabatic  pulsations  has 
been  applied  and  is  in  fairly  good  agreement  with  observation.  In  the  at- 
mosphere of  a  variable  star,  however,  amplitudes  become  large  and  non-linear 
terms  important.  This  is  confirmed  by  the  fact  that  the  observed  variations 
of  brightness  are  not  sinusoidal  functions  of  time  but  show  a  pronounced 
skewness.  It  was  therefore  decided  to  find  numerical  solutions  for  the  original 
non-linear  equations  of  hydrodynamics  but  with  the  following  simplifying  assumpt- 
ions: 

(i)  The  Cepheid  atmosphere  is  considered  to  be  plane -parallel. 

(ii)  The  atmosphere  is  originally  isothermal  and  in  hydrostatic 
equilibrium  under  constant  gravity. 

(iii)  All  changes  are  adiabatic. 

(iv)  The  influence  of  the  interior  of  the  star  is  approximated 

by  a  sinusoidal  motion  of  the  bottom  lajrer  of  the  atmosphere. 


1)  cf.  Lautman,  Don  A.,  Doctoral  Thesis,  Princeton  University,  1956.  To  be 

published  in  Astrophys.  J.,  1957* 

2)  cf.  Eddington,  A.  S.,  The  Internal  Constitution  of  Stars,  chapter  viii. 

Cambridge  University  Press,  1926. 
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The  formulation  of  boundary  conditions  for  the  top  of  the  atmosphere 
is  a  much  more  difficult  prol)lem.  An  isothermal  atmosphere  theoretically 
extends  to  infinity,  but  in  a  numerical  computation  it  clearly  must  be  cut  off 
at  a  certain  point.  The  upper  boundary  condition  should,  therefore,  express 
the  effect  of  the  neglected  part  of  the  atmosphere. 

It  is  essential  for  this  investigation  that  the  oscillation  of  the 
bottom  layer  generates  outgoing  traveling  waves  of  higher  frequency  which 
are  not  reflected  at  the  "surface"  of  the  atmosphere.  The  upper  boundary 
condition  was  constructed  to  fulfill  this  requirement. 

21.21  Basic  Equatioas  of  the  Problem 

We  introduce  the  following  symbols 5 


Y     velocity 
t      time 
g      gravity 


r    position      r^  A    same  as  r,  f     ,  p, 
o    density    j   O©  r       ^'^^  ^°^   hydrostatic 
p    pressure      P   J    equilibrium 


The  time  t  aa3  the  "iaitial"  petition  r^  are  used  9.a   ind«.pendent  variables. 
The  four  basic  equations  thsr.  raad; 


Equation  of  motion  — r-  =  "  '^  "  ~  '^~~^^  ^   ' 


Hydrostatic  equilibrium      ■— -  =  -  Opf^  (2) 

Equation  of  continuity       ^-    —      (■37]  '■^^ 

Adiabatic  equation  -^  =  \-%~\  ^   ' 


If  p/p  is  used  as  variable,  eq.  (1)  can  be  written  in  the  form 
At  T-o      o\1     V  pa.'    0  \  P«     / 


(5) 
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where  eq.  (2)  has  been  used  to  eliminate  dp  /dr  .  This  transformation  would 
simplify  the  computation  since  implicit  use  of  eq.(2)  eliminates  the  necessity 
to  compute  the  erponential  function  for  p  and  p  •  It  also  avoids  certain 
scaling  difficulties  on  a  fixed-point  computer. 

The  boundary  condition  at  the  bottom  is  simply 

V  =  AsincOt      .  (6) 

A  suitable  boundary  condition  for  the  top  can  be  derived  from  the  analytic 
solution  of  the  linearized  equations. 

21.22  Analytic  Solution  of  the  Linearized  Equations 

Let  us  put  u=  r-r  for  the  displacement  from  the  position  of  static 
equilibrium,  hence  v  =dr/ht   =3u)/<3t  .  From  equations  (3)  and  {k)   we  obtain 

We  now  use  the  assumption  that  the  oscillations  have  small  amplitudes,  thus 

This  brings  eq. (5)  into  the  linear  form 

where  c  =  po/^o  is  a  constant  for  ideal  gases.  Assuming  a  time  dependence  of 
the  form  ezp  (iQt),  the  solution  will  greatly  depend  on  the  frequency  D.  .   If 
a    is  very  high,  the  wave  length  is  very  short  and  the  gravitational  term  can 
be  neglected.  These  short  waves  are,  therefore,  traveling  with  the  speed  of 
sound,  viz.  -^i^ 

""-"of 

?or  a  next  approximation  let  us  put 

y  =  exp  (ikr  -iQt) 


p-p  I  «  p  and 

1*^  *^ol    *^o 


1 

P-R5 

t 

p 

r 

3l^ 
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hence  d  /^r^  =  Ik  and  c>/dt=s  —iSl  ,  Introducing  these  expressions  in  (7)  we  get 

or,   if  k  »  g/c^ 

Replacing    I'D.  and  ik  by  the  correspanelinc^  operators  we  get 

i!£_  «=,  ^ ^ iil  (8) 

brc  2d-       Co^'/i  dt     ^ 

which  was  used  as  upper  boundary  condition  for  most  cases.  As  a  variation,  a 
slightly  different  boundary  condition  was  used  to  represent  a  high  temperature 
layer  (corona)  at  the  top  of  the  atmosphere. 

21.23  Computational  Procedure 

As  long  as  the  linearized  equations  can  be  used,  the  harmonic  perturba- 
tion at  the  bottom  layer  will  generate  a  harmonic  movement  of  the  whole  at- 
mosphere. In  the  general,  i.e.  non-linear  case,  this  movement  will  atill  be 
periodic,  but  not  harmonic.  It  was  felt  that  the  periodicity  could  not  be 
used  efficiently  to  simplify  the  problem.  The  time  t  will,  therefore,  appear 
explicitely  as  an  independent  variable  together  with  r^. 

The  computation  was  started  from  the  hydrostatic  equilibrium,  i.e. 

r  a  r  and  v  =  0  with  an  oscillation  (cf .  eq.  6)  at  the  bottom.  The  integra- 

o 
tion  of  the  differential  equations  was  continued  until  the  motion  of  the  whole 

atmosphere  became  essentially  periodic. 

In  order  to  perform  the  numerical  computations,  the  (r^,  t)  plane  is 
divided  into  a  grid  of  constant  spacing  4  r  and  /^t.  To  obtain  the  highest 
accuracy  from  simple,  low-order  approximations  it  is  necessary  to  use  centered 
differences.  This  requires  that  the  variables  be  carried  on  staggered  points 
of  the  grid  according  to  the  following  list: 
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VARIABLES  SPACE  STEP  TIME  STEP 

y,  ^y/<^t  integer     ....    integer 

V  =  ^^/dt  integer  .    half -integer 

^\)/dr^   ,  p,^  -  .    .   half-integer    .   .   •    integer 

All  variables  are  then  replaced  by  dimensionleaa  variables  such  as  v  =  v/c 
and  z  =  tg/c  for  speed  and  time,  respectively,  and  all  differentials  are 
replaced  by  simple  central  differences. 

The  computational  procedure  is  as  follows;  Suppose  we  know  all  relevant 

variables,  including  ay/6t,   for  a  certain  time  t  and  that  our  time  step  equals 

At.     First  v(r  )  is  computed  for  all  r  's  from 
^  o  o 

except  at  the  boundaries  where  eq.(6)  and  (7)  are  used.  Similarly 


v(t-  4L)  ^  At.  4f  (t)  ===>  v(t  +  ^)   , 


y(t)  +  ^t.  v(t  +  ^^)  ==>  \j  {t  +  At) 

Finally;,  the  finite  difference  equivalents  of  eqs.  (3),  (k) ,   and  (1)  or  (5)  are 
used,  in  this  sequence,  to  compute  ^  ,  p  and  c>v/c)t  for  t  +A  t,   whereby  the  cycle 
is  completed. 

One  computational  detail  may  be  worth  mentioning!   In  eq.(i)-)  <? /pc   must 
be  raised  to  the  power  V"  ,  which  in  this  computation  was  taken  as  5/3.  Original- 
ly, an  iterative  method  with  constant  initial  "guess"  was  used  to  compute  the 
cube  root,  but  it  was  soon  discovered  that  this  took  almost  half  of  the  tetal 
computation  time.   It  was  replaced,  therefore,  by  a  simple  rational  approxima- 
tion, viz. 

5/3 
(Ah  +  B)'  z  +  C  ==>  z'"-' 

for  each  of  three  ranges  of  z,  yielding  an  accuracy  of  about  1<^. 

21.2l«-  Numerical  Stability 

When  approximating  a  differential  equation  by  a  difference  equation  we 
expect  certain  errors  in  the  solution.  One  type  of  error  (truncation  error) 

arises  from  the  neglect  of  higher  terms  in  the  series  expansions  which  replace 
the  derivatives.  The  truncation  error  can  usually  be  made  as  small  as  one  wants 
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by  either  using  a  higher  order  approximation  for  the  derivatives  or  by 
reducing  the  grid  spacing.  Another  type  of  error  is  due  to  the  fact  that 
the  difference  equation  may  not  be  stable  even  though  the  differential 
equation  is. 

In  order  to  avoid  this  kind  of  instability,  an  elaborate  analysis 
was  carried  out  for  whJch,  however,  our  system  of  equations  had  to  be  linear- 
ized. The  results  of  this  analysis  hold  strictly  only  for  the  linearized  system 
but  they  give  a  fair  indication  of  what  one  might  expect  from  the  non-linear 
equations. 

The  result  from  the  linearized  system  is  that  the  Courant-Priedrich 
1) 


condition 


>  4U,    =  velocity  of  sound 


At 

must  be  fulfilled  which,  for  a  given  space  step,  sets  an  upper  limit  for  the 
time  step.  But  because  the  atmosphere  will  be  compressed  during  part  of  the 
pulsation  cycle  and  the  velocity  of  sound,  therefore,  increased,  it  becomes 
necessary  to  make  At      smaller  than  A'r-^/w 

*      We  previously  mentioned  the  advantages  gained  in  simplicity  and  speed 
of  computation  if  eqo(l)  is  replaced  by  eq. (5) .   In  actual  computation,  however, 
it  is  found  that  the  numerical  solutions  of  this  modified  system  suffer  from 
a  very  strange  type  of  instability.  The  errors  increase  slowly  to  a  fairly 

constant  value  which,  though  small,  is  big  enough  to  destroy  the  periodicity 

2) 
of  the  solutions.  A  possible  explanation  is  given  in  Lautman's  Thesis   to 

which  we  refer  the  interested  reader. 

21.25  Results 

■      The  numerical  integration  was  carried  out  for  four  different  frequencies 
and  for  two  different  amplitudes  at  the  bottom,  while  the  height  of  the  atmosphere 
was  generally  taken  as  10  "scale -heights " ,  H  =  c   /g.   In  addition  to  these  eight 


1)  Courant  B,  Friedrichs  K.  and  Lewy  H.  ,  Math, Ann. 100;  page  132,(1928) 

2)  Lautman,  Don  Ao,  Doctoral  Thesis,  Princeton  University,  1956 
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runs,  several  experiments  with  varying  parameters  (height,  time  and/or  space 
steps)  were  carried  out,  and  also  one  involving  a  corona.  Each  run  took  about 
one  hour  on  the  machine. 

The  first  period  will  always  show  strong  transients  which  soon  settle 
down;  the  motion  is  essentially  periodic  by  the  sixth  period. 

The  results  show  that,  even  with  the  simple  model  which  we  have  chosen, 
it  is  possible  to  account  for  the  observed  skewness  of  the  velocity  curve  (as 
a  function  of  time)  and  for  certain  humps  on  the  descending  branch.  Both  of 
these  effects  are  very  strongly  dependent  on  the  boundary  condition,  however, 
and  in  further  work  on  this  subject  it  will  be  necessary  to  analyse  in  much 
more  detail  the  effect  of  the  upper  parts  of  the  atmosphere  on  the  character 
of  the  pulsation.  The  travelling  wave  components  which  lead  to  the  phase  lag 
are  seen  to  arise  quite  simply  from  the  standing  wave  fundamental  of  the  interior, 
but  the  magnitude  of  the  lag  is  slightly  less  than  that  observed.   It  should 
again  be  emphasized  that  the  simplifying  assumptions  which  have  been  made 
--  notably  the  adiabatic  approximation  and  the  plane -parallel  atmosphere  -- 
probably  influence  the  final  results  greatly. 


1) 
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21.30  MAGKETOHYDBODYITAMIC  THEORY  OF  SOIAR  SPICULES 
Originator  and  Coder;  Eeimar  Lust 

During  the  last  two  years,  calculations  concerning  this  problem 
have  been  carried  out  by  Dr.  Schliiter,  Dr-  Schwarzschild  and  Mr.  Lautman 
of  the  Princeton  University  Observatory.  They  have  calculated  pressure 
waves  which  expand  in  an  isothermal  atmosphere  from  the  base  of  the 
chromosphere  in  the  presence  of  a  vertical  magnetic  field.  An  essential 
simplification  for  this  problem  was  the  assumption  of  waves  with  small 
amplitudes  which  allowed  the  basic  differential  equations  to  be  linearized. 

The  aim  of  the  present  investigation  is  to  calculate  waves  of  arbit- 
rary amplitude.  Therefore,  the  non-linear  differential  equations  have  to 
be  integrated.   In  the  equations  of  motion  an  artificial  viscosity  term  was 
introduced  in  analogy  to  the  method  of  von  Neumann  and  Eichtmyer  'in  order  to 
allow  for  the  possibility  of  considering  the  transition  of  a  wave  with  finite 
amplitude  into  a  shock  wave. 

21.31  Basic  Equations 

The  basic  equations  for  this  problem  are  the  following? 

Momentum  equation: 

(1)  ^^  =  ]:*.  H    -  c^rccclC  P+q)  -  ^^  c^rv^\<^. 


Equation  of  continuity; 


ct 


(2) 

Energy  equation;  ^p  _  f   P    .    ,  3  "l  dg 

(3)  ^t  -  i  ^7^  ^^  '^_p  J  ,(t 

Ohm's  law  (with  (S'  =  oc   ): 

i^)  E  =  -  7^  H 


■'"^     For  detailed  description  of  the  astrophysical  problem  of  spicules  see 
Final  Report  on  Contract  No.  DA-36-03^-aRD-1330,  11,56  (195^)- 

2)     von  Neumann,  J.  and  Richtmyer,  R.  D. ,  A  method  for  the  numerical 

calculation  of  hydrodynamic  shocks.  J.  Appl.Phys.  21,  232-237  (1950). 
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Maxwe 11 ' s  equat  iona ; 

(5a)  VxE     =    -    9t 


(5b)  v-H    =  4tuJ  ('^    IIt'^^'I'^ 


where 


o   =  density,  p  =  pressure,  v  =  velocity,  ^  =   gravity  potential, 

E  =  electric  field,  H  =  magnetic  field,  j  =  current  density, 

>r  =  ratio  of  the  specific  heats,  t  =  time. 

q  is  an  artificial  viscosity  term  and  is  given  by  the  fcllcwing 

expression; 

[3  is  a  constant,  and  .^  x  is  the  interval  length  in  space  used  in 
the  numerical  computation. 

21.32  Assumptions 

The  same  assumptions  as  in  the  previous  investigation  have  been  made^ 

namely, 

a)  plane  parallel  chromosphere, 

b)  barometric  density  gradient    (  p  ---  e    ;  ^=  9-J>f>.=  Ccrv_  .) 

c)  small  vertical  extent   (  VCJ)  =  g  =  const.), 

d)  axial  symmetry  (H^  =  v^^  =0;  Eh  =  Jz  =  Js  =  °)  y 

e)  homogeneous  vertical  magnetic  field  of  value  H^, 

f)  all  changes  adiabatic,  apart  from  the  artificial  viscosity  term. 

21.33  Axially-Symmetrical  and  N  on -Dimensional  Equations 

With  the  assumption  of  axial  symmetry  all  independent  variables  should 
be  only  a  function  of  the  two  space  variables  s  (distance  from  the  vertical 
axis  of  symmetry)  and  z  (distance  from  the  horizontal  base  plane)  and  of  the 
time  t.  The  velocity  v  and  the  magnetic  f ield  H  are  given  by  the  horizontal 
components  v^  and  H^  respectively,  and  by  the  vertical  components  v^  and  H^- 
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The  following  non-dimensional  variables  are  introduced; 

Dq  =  undisturbed  density;  p  =  undisturbed  pressure. 

With  these  variables  one  gets  the  following  equatiorB  after  eliminat- 
ing  the  electric  field  E  and  the  current  density  j:  and  wr^Kn^j  f/p,  <^i  ■\<^^   ^.i°/M 

ax        dr  ^'-^,       ^   ^   "(^        or  cij  / 


(6a) 
(6b) 
(7) 


'''      t    --^t--^t^--[^P^^^-'^^](^^-|7-a^-)-^Pv 


(9a)   AS   _.„  _v^_S9v^^|2^2|^ 

I  dr   ^     ^^  ^ 


(10,   ^   =  .p(,,3-p(^.|^,|^),|^_.|_w,^^ 


21. S**-  Boundary  and  l-nitial  conditions; 

To  simplify  the  problem  it  is  assumed  that  the  wave  is  travelling  in 
a  "box".  Therefore,  the  normal  components  of  the  velocity  shall  vanish  at 
the  boundaries.  For  the  top  the  assumption  is  made  that  both  components  of 
the  velocity  are  zero  in  order  to  avoid  incoming  waves.  Besides,  boundary 
conditions  for  the  magnetic  field  have  to  be  introduced.  As  a  simplification 
it  is  assumed  that  the  magnetic  field  is  always  perpendicular  to  the  bottom, 
and  parallel  to  the  vertical  boundaries. 

A3  initial  condition  a  "pressure  hump"  is  assumed  which  should  start 
from  the  origin.  The  initial  values  for  the  6  variables  are  then  given  by 
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w  =  V  =  0 

p=^=  1  +  C.exp[-a*(At^^)] 

S  =  0 

Z  =  1. 

Here,  C  gives  the  amplitude  of  the  pressure  hump  and  oc^  ia  inverse 
proportional  to  the  decay  length. 

21.35  Difference  Equations 

To  convert  the  differential  equations  (6a)  -  (10)  into  difference 
equations,  the  space  differential  quotients  are  replaced  by  difference 
quotients  correct  to  the  second  order,  while  the  time  differential  quotients 
are  replaced  by  difference  quotients  correct  to  the  first  order;  r,  y,  and 
Z    are  introduced  as  independent  variables • 

T     =    -R  A  r  0<E<iP 

y    =    T>4y  0  <    Y  <   Tl 

The  differential  quotients  are  then  given  by 


'cT  |,^Y  "  2Arr''?*'.Y-^^i^-i,Y 


)  '     5^U^4^^^'--"-'-^ 


and  analogous  for  the  other  variables. 

21-36  Coding  and  Numerical  Computation8_ 

The  problem  could  not  have  been  treated  without  the  aid  of  the  new  large 
magnetic  drum.  With  this  appreciably  enlarged  storage  it  was  possible  to  in- 
crease the  number  of  points  in  the  grid,  which  was  chosen  to  be  19  x  30  points, 
(J5  =  ^9f   Xl'  ~  ^^^ '  ^^  ^^^   drum,  the  values  of  the  6  variables  w,  v,  ^  ,   p, 
S,  and  Z  for  each  point  from  the  preceding  time  step  T  were  stored.  Further, 
the  values  of  the  next  time  step  T  +  1  were  stored  at  those  storage  places  where 
the  values  of  the  time  step  T  -  1  had  been  before.  The  values  for  one  point 
were  stored  in  successive  places.  The  values  for  three  lines  of  the  grid,  namely 
for  R  -  1,  B,  and  E  +  1  are  needed  for  the  calculation  of  the  line  R  for  the  next 


21.36 

time  step.  Since  the  Williams  memory  was  not  large  enough  to  store  all  the 
values  of  three  lines,  each  line  was  split  into  three  parts  which  were  separate- 
ly shifted  from  the  drum  to  the  Williams  memory  and  back. 

The  whole  code  was  also  stored  on  the  drum  because  the  Williams  memory  was 
too  small.  The  code  for  the  six  difference  equations  was  divided  into  three 
parts;  part  one  for  the  calculation  of  the  values  at  the  bottom,  part  two 
for  the  main  field  including  the  two  vertical  boundaries,  and  part  three  for 
the  top.  Further  parts  of  the  code  contaiced  the  calculation  of  the  initial 
values  and  the  orders  for  punching  out  the  results. 

For  reasons  of  checking  every  time  step  was  calculated  twice.  The  cal- 
culations were  only  carried  on  if  both  results  were  in  ag;:eem.eiito   if  the 
results  differed,  the  machine  stopped  to  provide  information  for  the  operator 
to  the  effect  that  an  error  had  been  made.  Eestarting  caused  the  time  step 
to  be  repeated  again.  The  integration  of  one  time  sxep  for  all  the  points  of 
the  grid  took  about  two  minutes.  The  magnitude  of  a  "cime  interval  determined 
how  often  the  results  were  punched  out,  and  the  values  for  every  second 
point  of  the  grid  were  punched.  The  time  needed  for  punching  out  the  values 
for  one  time  step  was  about  k   1/2  minutes.  Extensive  restart  procedures 
were  incorporated  into  the  code  so  that,  in  principle,  no  more  than  about 
10  minutes  of  computing  time  would  be  lost  due  to  machine  errors. 

21.37  Results 

So  far,  for  two  sets  of  parameters  lUo  time  steps  have  been  calculated. 
In  the  first  case,  the  amplitude  of  the  pressure  hump  was  small  (c  =  0.01) 
and  resulted  in  waves  with  small  amplitudes.   In  the  second  case,  the  amplitude 
of  the  pressure  hump  was  equal  to  the  initial  pressure  (C  =  1)  and  led  to  waves 
with  much  larger  amplitudes.   In  both  cases,  the  strength  of  the  magnetic  field 
was  chosen  in  such  a  way  that  the  magnetic  pressure  was  equal  to  the  gas  pressure 
in  the  middle  between  the  bottom  and  the  top.  Though  according  to  the  computa- 
tions which  have  been  carried  out  until  now  the  wave  has  not  yet  reached  the 
top  of  the  grid,  one  can  see  the  influence  of  the  magnetic  field.  The  magnetic 
field  has  the  effect  that  the  waves  are  strongly  guided  into  the  vertical  direct- 
ion. 
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2lA0  AN  IMPROVED  SOLAE  MODEL  WITH  THE  CAEBON  CYCLE  INCLUDED 
Originator  and  Coder:  Ray  Weymann 

A  recent  model  of  the  solar  interior  by  Schwarzschild,  Howard,  and 
Harm  ^  in  vhich  all  energy  was  generated  by  the  proton-proton  reaction, 
appropriate  for  dwarf  stars,  indicated  that  the  central  temperature  was 
approaching  the  point  where  the  carbon-cycle  (which  provides  the  energy  for 
more  massive  stars)  would  contribute  significantly  to  the  total  amount  of 
energy  generated.   In  the  present  study,  the  effect  of  the  carbon  cycle  was 
taken  into  account.   In  addition,  an  improved  distribution  of  the  variable 
chemical  composition  of  the  sun  was  used. 

The  equations  involved  are  those  of  mechanical  and  thermal  equilibrium, 

with  two-point  boundary  conditions,  determining  two  eigenvalues.  These  are 

the  first-order  non-linear  ordinary  differential  equations,  whose  form  varies 

in  different  parts  of  the  star.   In  the  outer  portion,  the  envelope,  the 

2) 
equations  are  those  appropriate  to  convective  equilibrium  (cf.   and  Section 

21.5  of  this  report,  eq.  (1)),  while  in  the  central  region,  the  equation  for 
the  temperature  distribution  is  that  appropriate  for  radiative  equilibrium. 
Since  the  envelopes  had  been  previously  integrated,  it  was  now  necessary  to 
integrate  the  equations  for  the  central  region  and  fit  the  two  solutions 
together  at  the  point  at  which  convection  seta  in. 

The  entire  series  of  integrations,  including  computations  of  the  start- 
ing values  (which  must  be  done  by  series,  due  to  the  existence  of  singularities 
in  the  coefficients  at  the  center)  was  performed  on  the  computer,  using  the 
FLINT  code,  and  was  fitted  to  the  envelopes  by  hand.  The  results  for  this 
improved  solar  model  differ  only  slightly  from  those  obtained  previously. 


■'•^      Schwarzschild,  Howard,  and  Harm,  Afltrophys.  J.,  125,  233(1957) 
2)      Osterbrock,  Astrophys.  J.,  II8,  529(1953). 
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21.50  ON  THE  PULSATIONAL  STABILITY  OF  STAES  WITH  CCNVECTIYE  ENVELOPES 

Originator  and  Coder;  Irving  N.  Eabinowitz 

In  order  to  study  the  possibility  of  cepheid  pulsations  being  main- 
tained by  thermonuclear  energy  generation  in  the  deep  interior  of  the  stars, 
it  is  necessary^  first  to  construct  models  of  cepheid  variables,  and  second, 
to  solve  the  pulsation  equations  arising  from  a  perturbation  of  the  model. 
In  this  investigation,  three  models  were  constructed  which  differed  from 
those  studied  by  previous  investigators  in  that  they  had  convective  envelopes 
rather  than  radiative  ones.  Such  envelopes  are  described  by  the  equations 

These  are  the  equations  for  conservation  of  mass,  mechanical  equilibrium,  and 
convective  equilibrium,  respectively.   Integrations  of  these  equations  had 
been  previously  performed  by  Harm  and  Schwarzschild   on  the  IAS  computer, 
but  not  to  depths  great  enough  for  the  present  study.  It  was  therefore 
necessary  to  extend  the  integrations  for  the  three  envelopes  which  were  used. 

The  wave  equation  derived  from  linearizing  the  general  perturbation 
equation  was  broken  down  into  two  first -order  equations  for  the  radius  change 
and  the  density  change,  respectively: 

g\  ?'-^3r'   =0  ^  (2) 


|f'-V(y)p'  -Y{y)[2.k.^lc^^] 


r'  =  0 


where  y  =  logx  . 

The  quantity  CO  ^  is  an  eigenvalue  representing  the  period  of  pulsation, 
and  these  equations  have  well-behaved  solutions  for  only  discrete  values  of  OJ   . 


1)  cf.  Osterbrock,  Astrophys.  J.,  Il8,  529(1953)' 

2)  Harm  and  Schwarzschild,  Astrophys.  J.Suppl.  Series,  1,  319(1955) 
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21.51  Methods  of  Computation 

The  method  of  determination  of  the  eigenvalue  was  as  follows:  several 
integrations  of  (2)  for  various  values  of  OJ  were  done  on  the  machine  and 

plotted  by  hand.  From  the  graphs,  it  was  possible  to  choose  another  set  of 

2 
CO  values  more  closely  approximating  the  desired  value.  This  process  con- 
verges quite  rapidly,  since  by  interpolation  in  the  graphs  it  is  possible  to 
guess  at  the  eigenvalue  with  perhaps  one-and-a-half  decimal  place  accuracy 
in  the  differences,  so  that  no  more  than  six  runs  on  the  machine,  each  requir- 
ing about  an  hour,  were  necessary  to  determine  the  eigenvalue  to  nine  places. 
This  procedure  was  carried  out  separately  for  the  three  models  in  question. 

The  method  of  integration  used  was  an  extremely  simple  one,  namely 
Heun's  method,  in  which  the  value  of  the  function  at  the  forward  point  is 
determined  by  the  trapezoidal  rule,  and  is  then  improved  by  taking  the  average 
of  the  two  slopes  and  using  this  to  determine  a  new  va].ue  at  the  forward  point. 
This  is  a  primitive  type  of  Eunge-Kutta  method,  of  second -order  accuracy  and 
rather  high  speed.   It  may  be  remarked  that  the  small  accuracy  of  the  integra- 
tion method  does  not  conflict  with  the  necessity  of  finding  the  eigenvalue  to 

full  accuracy,  the  reason  being  that  it  is  the  shape  of  the  density- change 

2 
curve  that  counts,  rather  than  the  value  of  CO  .  Thus,  changing  the  truncation 

error  by  using  a  finer  grid,  or  a  more  accurate  integration  method  has  the 

2 
effect  of  changing  the  value  of  oj  to  which  the  well-behaved  solution  belongs, 

but  has  practically  no  effect  upon  the  solution  itself. 

All  computations  were  performed  using  the  floating-point  interpreter 

(section  12),  and  took  a  total  of  about  twenty  hours  on  the  machine.  A  single 

integration,  including  punching  out  a  table  of  the  results,  took  about  four 

minutes . 

21.52  Results 

It  was  found  that  even  for  models  with  convective  envelopes  the  energy 
generation  by  thermonuclear  processes  in  the  deep  interior  was  insufficient 
to  maintain  the  pulsations,  since  dissipative  processes  in  the  envelope  lose 
more  energy  than  is  generated  as  a  result  of  the  pulsations.  The  same  property 
has  been  previously  found  in  models  with  radiative  envelopes.  Since  the  deep 
interior  is  incapable  of  explaining  the  long-time  stability  of  cepheid  variables, 
it  now  becomes  necessary  to  study  the  mechanisms  of  energy  storage  and  dissipa- 
tion in  the  envelope. 
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22.  ATOMIC  AWD  NUCLEAE  PHYSICS 


22.10  TEDS  dROUKD  STATE  OF  THE  HELIUM  ATOM 
Originator;  T.  Kinoshita 

Coder:  Mrs.  S.  Bargmann 

Prof.  Toiciiiro  Kinoshita  haa  carried  out  an  elaborate  and  very  accurate 
re computation  of  the  ground  state  of  the  Helium  atom  '    in  order  to  match  the 
improved  accuracy  of  measurement  of  this  energy  level. 

It  is  veil  knovm  that  this  three-body  problem  (nucleus  +  2  electrons) 
cannot  be  solved  analytically*^  Very  close  upper  bounds  can  be  computed,  how- 
ever, by  help  of  the  EITZ  variation  method  if  a  sufficiently  large  set  of 
suitably  chosen  coordinate  functions  is  used.  Two  small  corrections  --  one 
for  the  mass  polarization  of  the  nucleus,  the  other  for  relativistic  terms  — 
may  be  estimated  with  sufficient  accuracy  after  the  unrelativiatic  three-body 
problem  has  been  solved. 

Prof.  Kinoshita  used  more  general,  and  obviously  more  suitable,  co- 
ordinate functions  than  his  predecessors  and  he  finally  increased  the  number 
of  coordinate  functions  to  39.  The  major  part  of  this  huge  numerical  work 
was  carried  out  on  the  AEC-UNIVAC  at  New  York  University.  Preliminary  computa- 
tions, with  up  to  10  functions,  have  been  done  on  our  machine  in  Princeton. 
We  wish  to  give  a  short  description  of  the  mathematical  problem  and  the  method 
used  for  its  solution. 

22.11  Basic  Equations 

The  original  quantum-mechanical  problem  consists   in  minimizing  the 
integral 


E  = 

under  the  normalization  condition 

2 


(  -y  (Hv)^^ 


^=^11/  cloj  =  1 


Here,  If  is  the  (real)  wave  function  of  the  Helium  atom  which  depends  on  the 
6  coordinates  of  the  two  electrons;  H  is  the  Hamiltonian  operator,  and  d.oJ 


1)  Kinoshita,  T.  Phys.  RgV.  105,  1^90  (1957) 


I 


22.11 

indicates  iHtegration  over  the  6  coordinates.  Tkree  angular  variables 
are  separated  so  that  i»  effect  only  three  independent  variables  remain. 

If  ^   is  erpreased  as  a  linear  conbination  of  trial  functions  Xi 
with  coefficients  u.  to  be  determined , 

a 

y  =  22u-J-,  ,  (1) 

the  two  integrals  E  and  N  are  turned  into  quadratic  forms   in  the  coeffi- 
cienls  u.  t 

e'  =  ^cX'UiL'M  =  (uAu) 


n'  =  ZI./^^j^^^-^j  =  i^,'\^^) 


(u  is  the  vector  with  components  u. ,  A and  B  are  symmntric  matrices  with 
elements  OC.  .  and  p..  respectively,  and  B  is  positive  definite). 

The  original  problem  is  then  replaced  by  that  of  minimizing  e' 
under  the  eiibsidiary  condition  k'  =  1,  which  In  turn  leads  to  the  eigen- 
value problem 

Au.  —  \Bu  (2) 

The  lowest  eigenvalue  A  is  the  minimrom  of  e',  and  for  an  appropriate 
choice  of  the  trial  functions  Xi    it  furnishes  an  approximate  value 
(in  fact,  an  upper  bouad)  for  the  minimum  of  E.  With  the  help  of  eq.(l) 
the  corresponding  eigenvector  furnishes  aa  approximate  wave  function. 

In  OTir  computations,  n  was  equal  to  10. 

We  determined  all  the  eigenvalues  of  (2) ,  not  merely  the  lowest 
one,  and  all  eigenvectors. 

22 . 12  Computational  Procedure 

Although  the  algorithm  used  here  has  been  described  in  the  Final 
Eeport  on  Contract  Wo.  Da-36-03*<-CED-1330,  December  195^>  we  shall  short- 
ly indicate  it  here  for  the  reader's  convenience.  The  main  steps  are; 

1)  Diagonalization  of  B,  i.e.  determination  of  an  orthogonal  matrix 

U  such  that 

D  =  U*BU 


I 
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is  a  diagonal  matrix  (witla  positive  diagonal  elements  Si  ).  Here,  U*  is  tke 
transpose  of  U. 

2)  Formation  of  D~  '  ,  tlie  diagonal  matrix  with  elements  J''^'^, 

-l/p 

3)  Formation  of  the  matrix  product  UD  '  . 

h)   Formation  of  the  symmetric  matrix 

W  =  (ud"^/2)*a(ud"^/^) 

5)  Diagonalization  of  W  in  the  form 

A  =  v*w 

where  V  is  orthogcaal,  and  A  is  a  diagonal  matrix  with  elements  A^  . 

6)  Formation  of  the  matrix  product 

-1/2 
Z  =  (UD  ^''^)Y 

Then  the  X  .  are  the  eigenvalues,  and  the  columns  z.  of  the  matrix 
Z  are  the  eigenvectors  of  the  problem  (2) 

To  prove  this,  we  first  note  that 

*      */   -l/2v*   /   -1/2,      *     A 
Z  AZ  =  V  (UD  '  )  A(UD  '  )  V  =  V  W  =  A 

Z*BZ  =  Y%-'-^^   U*BUD-^/Sr  =  V*D-^/^  D"^/^  =  I 

Setting  u  =  Zy,  and  premultiplying  both  sides   of   (2)   by  Z     we   obtain 
the  equivalent  eigenvalue  problem 

Ay  =  Ay 

in  terms  of  the  vector  y..  Since  A  is  diagonal,  its  eigenvalues  are  X^, 
the  corresponding  eigenvectors  e.  (e.  has  the  i-th  component  1,  and  all 
other  components  0).  Thus  the  eigenvectors  of  (2)  are  Ze^  =  z^,  q.e.d. 

The  algorithm  described  here  requires  two  diagonalizations  which  are 
carried  out  according  to  the  Jacobi  method  (see  Final  Report  on  Contract  No. 
DA-36-03U-aRD-1023,p.II-50,  where  also  the  checking  procedures  are  described.). 
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We  succeeded  in  somewhat  improving  the  accuracy  obtained  before  by- 
increasing  the  number  of  rotations  and  by  testing  explicitly  for  the  best 
scaling  factor  at  intermediate  stages  of  the  computation,  readjusting  our 
numbers  accordingly. 
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22.20  BELATION  BETWEEN  THE  VIBRATION  FREQUENCIES  OF  A  CRYSTAL  AND  THE 
SCATTERING  OF  SLOW  NEUTRONS 

Originator:  G.  L.  Squires 

Analyst  and  Coder:  Hans  ,T,  Maehly 

*) 
Under  this  title  Dr-  G.  L.  Squires,  has  published  a  paper  '  reporting 

and  discussing  the  results  of  computations  which  were,  according  to  his 

specifications,  carried  out  on  the  Electronic  Computer  at  the  Institute  for 

Advanced  Study.  We  wish  to  supplement  this  paper  by  a  purely  mathematical 

description  of  the  problem  aad  an  outline  of  a  few  interesting  aspects  of 

eur  code. 

22.21  The  Computational  Problem 

The  entire  problem  can  readily  be  divided  into  3  parts. 

(i)  Computation  of  the  Frequencies 

Given  5  constants  a,  b,  c,  e,  f,  compute  the  eigenvalues  A  and  "frequencies' 

M  =  J  A   for  which  the  3x3  determinant 
/^m   ^   m 

I  T(a)  -y?I  I  =  0 

will  vanish,  where  I  is  the  unit  matrix  and 

T..  =  a  +  2b  -  aC.C,  -  bC,(C^  +  C,)  +  eS^  +  f(S.  +  S,) 

S.  =  sincx.  1    i  =  1(1)  3 


C.  =  coscx.   I  CK--  =  O(^)  ^-<5' 


*)  Phys.  Rev.  103,  30^-312  (1956) 
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By  making  full  use  of  the  various  symmetries,  the  number  of  grid  point* 
{Oi^}  oC^yOC-.)    ^^^  ^e  reduced  considerably,  namely  from  1000  to  I52  for  (f=   l8° 
and  from  27.000  to  2792  for  d'=  6°.  These  are  the  values  of  cT  for  which  the 
computation  was  actually  carried  out. 

(ii)  Computation  of  the  Frequency  Distribution  Function 

After  arranging  the  frequencies  Ufl,(0C,  ,  CXp^  OC,)  according  to  their  size^ 
determine  the  number  W.  of  such  frequencies  lying  between  u     and  n.    +AJU<     for 
some  given  and  finite  value  of  4;\  and  aa  a  function  of  la  .   --This  was  carried 
out  for  d"=6  ,  dividing  the  total  range  of  u.   into  approximately  250  equal 
intervals;  u     =  0   (Alk)    U    ^    ,      u         ~  250. 4M-  •  As  full  use  was  made  in  (i) 
of  the  inherent  symmetries  of  the  problem,  every  frequency  must  be  multiplied 
by  an  integer,  which  is  k8   inside  the  reduced  volume  and  correspondingly  smaller 
for  the  faces,  the  various  edges  and  vertices  .  Finally,  to  reduce  the  statistical 
fluctuations  arising  from  the  finite  size  of  (f   and  ^aa   ,   a  smoothing  procedure 
was  applied  to  the  "curve"  N.  =  N(yiA.)  (cf.  2.'d.2h) , 

(iii)  Averages  of  Various  Functions  over  N(  M^  ) 

After  the  distribution  function  N(  l/U  )  has  been  compated,  the  averages  of 
the  following  functions  have  to  be  found. 

A,  =  J-  ccf  h  2^  ] 


for  15  different 
values   of  CX^t 

V  = 

-2 

(i). 

2 

Since  Aj^   is  finite,  all  integrals  must  be  replaced  by  surna 


where 


N  =^n^  =  (  ^/<f  ) 


and  F  stands  for  Ay  ,  By  ,  or  C,  ,  respectively. 

Both  the  Ay's  and  the  first  C,  's  are  ini 
integrals  still  converge  as  n(M-)  tends  to  zero,  but  our  sums  are  becoming  very 


Both  the  Ay's  and  the  first  C,  's  are  infinitely  large  for  /M-=  0.  The 
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sensitive  to  and  greatly  affected  by  the  statistical  errors  of  n.  for  small 
i.  A  special  smoothing  procedure  had  to  he  developed  to  overcome  this  dif- 
ficulty  (see  22.2i<-)  . 

22.22  Square  Root  Subroutines 

Our  code  required  the  computalbion  of  very  many  square  roots.   It  was 
largely  due  to  the  development  of  fast  square  root  subroutines  that  the  run- 
ning time  for  this  problem  could  be  kept  within  quite  reasonable  limits.  The 
total  computing  time,  not  counting  in-  and  output,  for  oie  frsquency  distribu- 
tion function  ia  just  9  minutes.  During  this  time  2792  matrices  are  computed, 
their  eigenvalues  determined,  and  the  square  roots  of  these  eigenvalues  are 
taken  anddi3tr5.butsd  among -250  intervals  according  to  size.  The  total  number 
of  square  roots  computed  during  these  9  minutes  is  about  ^0,O0Q<, 

About  i<-0,000  of  them  had  to  be  evaluated  for  the  diagoaalisation  of  the 
matrix  T  by  the  Jacobi  method.  For  every  two-dimensional  rotation^ tan2cp  is 
first  obtained  from  which  coa2a)  ,  sin2ai  ,  costp  and  3in<^  must  then  be  computed. 
This  requires  the  computation  of  ^  1  +  (tan2c£>)    and  \j  1  +  cos2i5»   where 
|2(i>|<  ""-/),•  The  radicand  is  thus  known  to  be  in  a  quite  restricted  interval  in 
which  a  reasonably  good  linear  approximation  may  be  used  to  reduce  the  number  of 
iterations  to  only  two. ' 

We  shall  show  below  how  a  best  linear  approximation  for  |x  can  be  found 
for  a  <  X  <  b. 

Another  square  root  routine  was  written  for  part  (iii)  of  our  problem 
where  \lA,'      was  needed  for  equidistant  points  ^=A/U.m^  m  =  1,2,...,  ^njoy' 
If  m  is  not  too  small,  rm'will  be  a  reasonably  good  estimate  for  \l  m+1",  but 
a  much  better  one  can  be  easily  obtained  as  shown  in  (ii)  below.  As  only  little 


*)  The  same  technique  can  be  used  for  the  computation  of  square  roots  of  normal- 
ized floating  point  numbers,  or  of  any  number  after  sensing  and  counting  by 
how  many  places  the  radicand  may  be  shifted  to  the  left  without  overflow. 
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accuracy  waa  required  for  the  amalleat  values  ofyu  a  suigle  Iteration,  usiag 
the  formula  shown  below,  was  sufficient  in  our  case. 

(i)  Linear  Approximation  for  ^,     a  <  x  ^  b. 

Most  square  root  subroutines  are  based  on  the  iteration  formula 

-  (  w.  +  -  )   =>  w.  T 
2  ^  1   w.       i+1 

1 

It  is  well  known  that  w  — >  fp  for  i  — >  oo  and,  in  particular,  if 

w.  ={T^(l+£.),  £.«  1 
then  w^^^  a;  >r?  (1+  2  ^  ) 

IB  any  case,  w  (and  hence  w  ,  w,  ,  ...)  is  always  too  big  and  the  error  depends 
omly  on  the  relative  accuracy  of  w  . 

Let  us  assume  that  we  knew  lower  and  upper  beuads  for  x: 
a  <  X  <  b 

Then  the  best  choice  for  a  constant  initial  value  obviously  is  w^  =  (ab) 
The  error  will  have  its  maximum  at  the  ends  of  the  interval  where 

The  next  iteration  yields  for  these  points 


Nfa"      Pb 


2 


'/r  ^(»^'1 


IB  order  to  get  best  relative  accuracy  over  the  whole  interval  we  divide  by  the 

square  root  of  this  expression  viz. 

lAl ")  1/2 


e  =  (|[(V)^^(»^^1}^ 
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We  thus  have 

"2 

1 
~     29 

where 

^0  = 

1 
2© 

(  V.    +  -    )       =C+C, 
\  *   V,  y       o    1 


X 


■  (^"l'^'^  =l=IT  <^^'"'^'' 

The  naximum  error  ratio  is  Q   for  a  <  x  <  b.   It  can  easily  be  ahown  that 
this  is  the  beat  linear  approximation  for  ^^   ia  the  given  interval.   In 
our  cede  this  was  used  for  b/a  =  2,  but  the  method  pays  off  evea  better  for 
wider  intervals  as  the  following  table  shews. 

b/a  Q 

10  2.2i+722l 

10^  2.986556 

10^  3.9783'n 

10^  5-303391 

10  7.07l'*21 

(ii)   Algorithm  for  a  Table   of  fP 

This  algorithm  is  based  on  the  approximation  formula 


b/a 

9 

2 

1.007l*-98 

k 

1.02988!^ 

-10 

1.081809 

100 

1.318807 

LOOO 

1.703121 

\  1-^     l.'^/2     ^    2-' 


7^ 

which  holds  for  u  «  1  with  an  error  of  ^  ~  u^/lj-.  Substituting  l/(2ia+l)  for  u 
we  obtain 


with  a  relative  error  of 


A    z  h 


for  m  »  1 


a  +  ^/2 
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Numerical  example: 

vflo^^  jZ9 '(i  +    9J73-)    =  7.0710660 

exact  value       =  7.07IO678 

absolute  error      =  .OOOOOI8 

relative  error      =  .00000026 
A  standard  iteration  vould  bring  this  relative  error  down  to  one  half  ita 
square,  which  is  less  than  J+.IO    for  our  example. 

22.23  Checking  Procedures 

Particular  care  has  been  taken  to  detect  and  trace  errors  due  to  false 
coding  or  machine  trouble.  We  wish  to  outline  some  principles  without  going 
into  toorauch  detail. 

(i)  No  attempt  was  made  to  minimize  the  number  of  locations  used  for 
temporary  storage;  on  the  contrary,  as  many  different  locations  as  feasible 
were  used  for  the  various  variables  so  that  a  maximum  of  information  about 
the  immediate  past  was  available  at  any  given  moment. 

(ii)  A  subroutine  was  written  and  incorporated  which  would  reconvert  to 
decimal  and  punch  out  on  cards  the  whole  field  of  temporary  storage  ("Post  Mortem 
Boutine") . 

(iii)  Transfer  to  this  subroutine  was  ms-de  (a)  manually  during  debugging 
after  each  "portion"  of  newly  tested  code  (b)  automatically  if  one  of  the  math- 
ematical checks,  incorporated  in  the  code,  would  fail. 

The  most  important  mathematical  check  consisted  of  squaring  each  frequency 

ix  just  after  being  computed  and  to  test  the  vanishing  of  the  determinant  |  T  -/^  1 1  • 

-9 
If  the  value  was  found  to  be  bigger  than  10  ^ ,   transfer  to  the  Post  Mortem  Routine 

was  made.  This  happened  just  once,  during  production,  and  it  was  possible  to  trace 

the  error  to  a  multiplication,  though  some  hundred  operations  had  been  carried  out 

by  the  machine  before  the  error  was  detected  by  the  determinant  check.  Both  factors 

together  with  the  faulty  result  could  be  found,  thus  providing  the  maintenance 

engineer  with  accurate  information. 
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In  addition  to  these  automatic  checks  several  mathematical  checks  were  carried 
out  manually,  to  test  the  code  and  the  accuracy  of  the  various  approximation 
procedures.  For  example,  the  third  part  (see  22.21-(iii))  was  run  for  a  "distri- 
bution" n.  =  1  for  0  <  1  .<  .  <  a  and  n.  =  0  beyond  >u.  =  a.  The  values  of  <  F  > 
thus  obtained  were  compared  with  the  integrals  computed  analytically. 

22. 2U  Smoothing  and  Graphing 

Due  to  the  finite  number  of  points  for  which  the  eigenfrequencles  were 
computed  the  statistical  fluctuations  of  the  n.  were  rather  distuxblng  because 
they  concealed  the  true  character  of  the  distribution  f-inction.  Therefore,  a 
subroutine  was  written  which  would  replace  each  n  by 

tv  =  r  (n.  ,  +  2n.  +  n.  ^) 
'H    ¥  ^  1-1     1    i+l"* 

This  "smoothing"  could  be  repeated  an  arbitrary  number  of  times,  but  after  several 
applications,  and  before  the  "noiae"  had  virtually  disappeared,  the  characteristic 
and  desired  "true"  discontinuities  of  the  curve  were  lost.  The  easiest  and  most 
efficient  way  to  find  the  best  compromise  was  to  use  the  graphing  unit  which  had 
just  been  completed  for  our  machine.  The  n.,  n.  n^ ,  . = .  were  displayed  on  the 
screen  of  a  cathode-ra5>-  tube  by  help  of  an  auxiliary  code  and  the  optimum  amount 
of  smoothing  could  be  determined  by  Dr.  Squires  himself  by  visual  Inspection. 

A  much  more  serious  smoothing  problem  arose  in  part  (iii)  of  the  code.   It 
can  be  sLoT<m  from  the  power  series  development  of  the  matrix  T  that  the  power 
series  for  n(yU.)  ,  at  ^=  0,  can  contain  only  even  powers  of  yU,  , 

2      h 
n(^)  =  C^^  +  C^^     +  .... 

_2 
On  the  other  hand,  the  Ay  and  C _^   grow  like  ^      on  approaching  ^=  0.  Hence, 

the  integrand  theoretically  tends  to  a  constant,  but  numerically  It  does  not 

because  of  the  "noise"  connected  with  the  n..  This  noise  is  proportional  to -fa^, 

hence  toyu.  ;   when  multiplied  by  yU  , for  C_2,it  will  increase  like  yU^  for  small 

values  of  >Lt  . 


Some  of  these  difficulties  could  have  "been  overcome  hy  replacing 
the  n.  by  n(/A  )  computed  from  the  power  aeries.  The  amount  of  work 
--  coding  and  machine  tine  --  to  compute  the  coefficients  Cp^  o^  't^® 
power  series  is  tolerable  for  k  =  0  but  grows  rapidly  for  positive  k. 
The  second  difficulty  arisea  at  the  joining  point  of  the  analytic  curve 
n(*x)  with  the  original  set  of  points.  Let  ua  examine  it  more  closely. 

We  recall  that  the  eigenfrequencies  u  have  been  computed  for  the 
meshpoints  only.  Each  meahpoint  ''represents"  a  certain  volume,  in  gener- 
al a  cube.  Each  of  the  three  eigenfrequencies  of  the  matrix  will  be, 
for  all  points  of  the  cube,  pretty  close  to  those  for  ita  center,  i.e. 
the  meshpoint,  but  not  necessarily  close  enough  to  be  in  the  aame  in- 
terval {u     ±  a'^)   ^^  ^^   corresponding  eigenfrequency  for  the  mesh- 
point. This  is  -'ihe   reason  for  the  "noise"  in  our  curve  n^. 

Let  ua  now  look  at  our  intergrala  which  are  represented  by  sums 
(see  22.21. iii).  The  straight  intergral 


\  r^f^ld/^  =   5;j'ZL*^^  —  I 


is  not  affected  by  the  f luctuationa ,  as  the  total  count  '^loes  not  depend 

on  the  esacx  location  where  the  eigenvalues  are  listed.     If,  however, 

2 
n(  ix)    is  !Bultiplied  b-w  aome  function  as  e.g.   f(yCA-)   =AA.  ,  then 


\  rvfyuUXhi^  "^  j:^  ^  "^v^' 


will  hold  only  approxissately.  But  even  in  this  case,  the  effect  of  having 
a  number  n  of  counts  erroneously  shifted  to  the  next  interval  will  cause  a 
relatively  small  error,  viz. 

Let  us  go  back  now  to  our  idea  of  replacing  the  n^  by  some  analytic  ex- 
pression n(^^)  up  to,  say,  /a=/\  and  let  us  assume  again  that  a  number 
n  of  counts  should  be  in  the  interval(/^±^<^but  are,  by  error,  in  the 
next  one.  In  this  case  the  error  will  be 
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rtf/Au-t-c^)  -^  ryxk 

since  ru    has  been  replaced  by  the  fixed  value  n(yU..  )  computed  from 
the  power  series.  This  error  is  much  bigger  than  the  first  one  since 

Therefore,  it  was  decided  to  join  the  analytic  function  n(/U^)  to 
the  "statistical  points"  n.  not  in  one  point  but  gradually  by  blending; 
the  new  values  r.  were  thus  computed  as 

rii  =  a-  ^(/aJ  +  ( i-(^ )  ru 

where  the  "weight  function"  g(/><-)  ia  defined  by 
and  ^      q_      —       i^ 

C  was  computed  directly  by  help  of  a  modified  code  replacing  each  matrix 
eleinent  by  the  first  term  of  its  respective  power  series  in  o( -if   ^o'  '^ V 

c^,  however,  is  a  "mean  value"  computed  from  the  original  n. (0  <  i  <  k)  by 

2'  Scune  ■  1   -   - 

the  least  square  method  using  the, weight  fu/action  gCu)  as  above. 

A  more  detailed  anfll3reis  has  shown  that  this  whole  procedure  should 
reduce  the  error  for  the  first  part  ;  0  < /(A  <  Aa,  (actually  about  a  quarter 
of  the  whole  interval)  to  roughly  the  same  amount  as  for  the  rest  of  the 
interval. 


The  development  of  this  smoothing  process,  together  with  the  special 
for  computing  C  ,  was  rather  time  consuming,  ye 
to  obtain  the  final  answer  with  sufficient  accuracy. 


code  for  computing  C  ,  was  rather  time  consuming,  yet  essential  in  order 
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22.30     NDMEEICAL  CALCULATIONS  OF  THE  ANOTLAE  DISTRIBUTIONS  FOR  THE 
DEUTEEON-PROTON  AND  SIMILAR  REACTIONS. 


Originator 

AnaljTBist 

Coders 


William  Tobocman 

Hans  J.  Maehly 

Sonja  Bargmann  and  Patricia  Eberlein 


When  a   light  nucleus    (e.g.   a  deuteron)    hits   a  heavier   one  with  suf- 
ficient energy  to  overcome  the  repellent  Coulomb   force,    i.e.   to  enter   into 
the  range  of  nuclear   interactions,   one  part  of  the   light  nucleus    (e.g.   a 
neutron)   may  be  captured  by  the  target  nucleus,  while  the  rest   (the  proton) 
will  leave,   accelerated  by  the  repellent  force,  and  at  an  angle  9  with  res- 
pect to  the  direction  of  the   incoming  particle.      In  the  experiment,  a  virtual- 
ly unidirectional  and  mono-ea^rgetlc  beam  of  light  nuclfei  is  aimed  at  a 
target  of  heavier  one3  and  the  density  (5(  0  )    of  outcoming  residual  particles 
(protons)    is   observed. 

This    "differential  cross   section"  <6i  Q  )   could  be  computed  if  the 
nature  and   laws   of  nuclear  forces  were  known;   as  they  are  not,  trial  and 
error  methods  are  used  to  confirm,    improve  or  reject   initial  guesses  for 
these  forces.     The  numerical  computation  nec«ssary  to  determine  (i{  9  ) , 
even  according  to  a  sircplo  "g'.ieaa"   of  the  forces  are,  however,  very  tedious 
in  most  cases.     Rather  rough  approximations  have  been  used  to  reduce  the 
work  to  a  tolerable  amount,   and  error  eatiraatea  are  so  difficult  that  the 
conclusions  drawn  from  such  rough  computations  seem  doubtful  in  some  cases. 

To  overcome  these  difficulties,  Dr«   Tobocman  made  the  following 
suggestions      : 

(i)   Not  to  start  with  the  nuclear  forces  themselves,  but  with  their 
effects    (phase  shifts)   at  a  certain  distance  R  from  the  center 
of  the  target  nucleus,  where  R  is  approximately  equal  to  its 
radius. 


■"■^  W.   Tobocman  and  M.  H.   Kaloa,  Numerical  Calculation  of   (d,p)  Angular 
Distributions.  Phys.  Rev.   97  (1955),   l,p. 132-136. 
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(ii)  To  prepare  a  moat  general  and  very  flexible  code  vhich  can 
be  used,  after  its  completion,  for  a  wide  range  of  initial 
data  (describing  the  experiment)  and  of  theoretical  aaaumpt- 
ion». 
Work  for  this  code  began  in  summer  1955  vith  an  extensive  study 
of  the  various  mathematical  problems ,  Programming  and  coding  waa  done 
through  the  year  1956.   The  code  is  now  in  the  debugging  and  testing 
phase.  It  is  expected  to  be  in  full  operation  in  summer  1957= 

About  SQfjo   of  this  work  was  supported  by  Contract  No.  Wonr-1358-(03) , 
the  rest  by  the  Army  Contract  for  which  this  volume  constitutes  the  Final 
Report.  We  shall  therefore  prepare  a  Technical  Report  to  the  Office  of 
Naval  Research  at  the  completion  of  the  problem. 
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22.i+0     DISTRIBUTION  OF  EIGENVALUES  OF  BORDEEED  MATRICES  WITH  INFINITE 
DIMENSIONS 

Originator:     Eugena  P.  Wlgner 

The  problem  to  determina  the  statistical  properties   of  the  character- 
istic values   of  infinite  bordered  matrices  arose  from  the   consideration 
of  the  properties   of  the  wave  functions   of  quantum-mechanical  systems  which 
are  assumed  to  be  so  complex  that  statistical  considerations   can  be  applied 
to  them.    ■'      The  mathematical  problem  together  with  an  integral  equation  and 
asjrmptotic  formulae  for  the  distribution  of  the  eigenvalues  and  with  a  re- 
currence  formula  for  their  moments  are  given  in  a  paper  by  WIGNER   ' ,  to 
which  we  would   like  to  refer  the  reader   interested   in  Vopt  details   of  the 
following  paragraphs . 

Our  task  was  to  find  a  method  and  to  write  and  run  a  code  for  the  nu- 
merical computation  of  the  distribution  of  the  eigenvalues  for  various  values 
of  the  parameter  q   (cf  eqs.    (1)   and    (!)•)  below).     The  results   so  far  obtained 
are  not  fully  satisfactory,  but  other  urgent  work  prevented  ua   from  improving 
them. 

22.^1       Two  Possible  Approaches 

(i)   Method  of  the   Integral  Equation. 

WIGNEE    (log.cit.)   derived  the  following  non-lin&ar   integral  equation 
for  the  distribution    o(x):  ^ 

L         -to-'  X  —  2.  -coy 


where 


1)  J.  A.  Lane,  E.   G.  Thomas  and  E.   P.  Wigner,   Giant  Resonance  Interpreta- 
tion of  the  Nucleon-Nucleus   Interaction,  Phys.Rev.   98   (1955),   pp. 693-701  • 

2)  E.   P.  Wigner,    "Characteristic  Vectors   of  Bordered  Matrices  with  Infinite 
Dimensions,"  Annals   of  Mathematics  62   (1955),  ?•   5^8' 
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Not  much,    if  anything,   is  known  about  thie  kind  of  integral  equation. 
But   it  was  hoped  that  the  follcwing  iterative  procedure  might  yield  an 
adequate  approximation  to  the  distribution  function. 

We  start  with  a  trial  function  p  (x)  which  may  be  constructed  by 
interpolating  very  roughly  between  the  asymptotic  solutions  for  very 
small  and  very  large  values   of  q,  respectively,  these  asymptotic  solu- 
tions being  known  frsmWigner's  paper.     We  then  compute  E   (x)    from  (2) 
and    9i^-^^    from  (1),   substituting    D     and  R     for  o  and  E  on  the  right- 
hand  aide,  and  so  on.     The  integrals  are  of  course  replaced  by  sums 
(e.g.   using  Simpson's  rule)   and  the   limits   in  (1)  by  some  finite  values 
±X.     From  the  known  asymptotic  behavior  of    P (x) ,  X  can  easily  be  deter- 
mined so  that  p  (x)   and  E(x)   are  numerically  negligeable  for  x  >  X. 

The  results   of  this  approach  were  unsatisfactory.      It  was  believed 
that  this  was   due  to  some  basic  difficulty  in  the  procedure.     Eencc  an- 
other method,  described  below,  was  devalopped;  but  it  should  be  said 
that  meanwhile  errors  have  baen  found  in  the  firat  code  so  that,  after 
all,  the   iteration  method  described  above  may  work  if  properly  coded. 

(ii)   Method   of  the  Moinentfl 

The   second  method  doss  uot  use  the  integral  equations    (1) 
and    (2)  but  a  recurrence  formrala  foi*  the  moments   of  the  di*itribution  p(x) 
which  is  also  found  in  Wigner's  paper   (loc.   cit.): 


All  odd  moments  vanish,  as  P (x)    is  an  even  function.      It  is   easy  to 
compute  yu     for,   say,    V  =  O....5O  from  (h) ,  particularly  since  all  con- 
tributions are  positive.     However,  the  solution  of  the  moment  problem, 
i.e.   computing  p(x)   from  the  ytu     ,    is,   numerically,  a  very  unstable 
procediore.     We  finally  succeeded   in  getting  fair  results   for  several  values 
of  q,  but  the  method  would  fail  for  both  very  large  and  very  small  values 


(3) 


(1^) 
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of  this  parameter.  Fortunately  the  known  asymptotic  solutions  just  fill 
thia  gap. 

2c;. ^2  Numerical  Solution  of  the  Moment  Problem 

Our  general  approach  is  to  try  to  fit  to  the  unknown  function  a 
step  function  having  as  many  moments  as  possible  correct.  We  use  step 
functions  with  a  fized  length  of  step  (denoted  by  h)  and  proceed  in  the 
following  fashion:  for  our  first  approximation  we  take  a  step  function 
with  one  step  such  that  its  integral  is  equal  to  the  integral  of  the  un- 
known function.  Next  we  take  a  function  having  three  sttps  such  that  its 
integral  and  its  second  moment  are  correct.  For  our  n-th  approximation 
we  take  a  sjmimetric  step  function  hs-ving  2n+l  steps  sach  that  its  first 
2n  moments  and  its  integral  are  correct.  Thus,  at  each  step,  we  are  solv- 
ing n  linear  equations  in  n  unknowns.   If  this  sequence  of  functions  con- 
verged, the  limit  would  be  a  function  with  the  correct  moments,  and  there- 
fore the  unknown  function.  However,  in  practice  it  does  net  converge,  and 
therefore  we  will  stop  the  process  at  some  point  where  we  believe  to  have 
obtained  the  best  approximation  to  the  function.  To  det^raii*  this  point, 
we  use  the  following  rather  hsuristic  argu.ment. 

With  our  n-th  step  function  we  are  matching  the  fiir-Jt  2n  moments 
of  the  unknown  function  with  a  function  of  extent  approxiiau.tely  from 
-2nh  to  +2nh,  constant  over  intervals  of  length  h.  The  extent  of  our 
step  functions,  therefore,  goes  up  as  2n.  With  these  functions  we  try 
to  approximate  the  unknown  in  that  portion  of  its  domain  which  gives  a 
significant  contribution  to  the  first  2n  moments.  This  region  does  not 
in  general  go  up  linearly  with  n.   If  this  region  of  significance  is  much 
smaller  than  the  extent  of  the  step  function,  it  is  clear  that  in  effect 
we  have  not  enough  parameters  in  the  region  of  significance  to  graph  the 
function  from  2n  moments.  Thia  is  why  the  infinite  sequence  of  step  func- 
tions diverges,  for  we  have  good  reason  to  believe  that  the  region  of  sig- 
nificance of  our  unknown  function  increases  at  a  less  than  linear  rate. 

Hence  our  approach  is  to  try  and  find  a  step  function  with  2n+l  steps 
such  that  its  extent  is  approximately  the  same  as  the  region  of  significance 
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of  the  unknown  function  for  its  first  2n  moments.  At  the  same  time  ve 
want  to  make  n  as  large  as  is  feasible  and  h,  the  length  of  the  steps, 
small  enough  to  graph  the  function  accurately. 

We  start  by  choosing  a,  particular  h,  and  proceed  to  apply  the  method 
of  approximation  described  above.  But  instead  of  using  simple  step  func- 
tions s  defined  by 


s^  =  1  tf(.n-\}<l<(-n.^)   or(n-  |)<^<fn+i) 


Srv  =  0  otherwise, 

we  take   linear  combinations   of  these  with  the  property  that  the  first  2n-l 
moments  are  zero,  while  the  (2n)th  moment  is  not.     This  puts  the  system  of 
equations -that  we  solve  at  each  step-into  triangular  form.     A3   long  as  the 
region  of  significance   of  the  unknown  function  is   larger  than  the  exigent 
of  the  approximation,  the  correction  to  the    (2n)th  moment  will  be  positive, 
that   is  the    (2n)th  moment  of  the  approximation  will  be  smaller  than  that  of 
the  unknown  function.     As   soon  as  this   is  no  longer  the  case,  we  stop  and 
take  the  approximation  at  that  step  as   our  best  fit  for  our  particular  choice 
of  the  step  len^h  h. 

For  each  value   of  q  we  have  used  several  values   of  h,  trying  to  obtain 
as  high  a  density  of  values  as  possible   in  at  least  the   intermediate  range 
of  X.     The  results  were  especially  good   in  the  range  of  q  between  .1  and  1.0 
Since  the  value   of  h  could  not  be  taken  too  small  (there  seemed  to  be  no 
reasonable  behavior  at  all  for  very  small  h)  we  could  not  obtain  points  close 
to  zero.     However,    in  the  range  mentioned  the  points   obtained  seemed  suffi- 
cient to  give  a  good  picture  of  the  function. 

For  very  small  q  <  .1,  the  behavior  around  zero  is  much  more  active, 
the  functicm  acquires  a  steeper  and  steeper  peak  which  is   impossible  for  us 
to  graph  accurately  with  this  method.     For  q  >  1,   our  technique  works  poorly 
and  we   got  fewer  and  less  reliable,  as  well  as   less   likely  points. 
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22.50     MOLECULAB   IFTEGRALS 

Originator;  R.  C.  Sahni 
Analyst  and  Coder:  J.  W.  Cooley 

Many  numerical  computations  have  been  carried  out  to  determine  the 
energy  levels  of  atoms  from  the  SCHROEDINGER  equation  and  the  relevant 
atomic  constants,  and  the  comparison  of  these  accurately  computed  theoretic- 
al values  with  the  observed  atomic  spectra  has  greatly  advanced  our  know- 
ledge of  the  laws  governing  the  electron  shellj  through  tne  analysis  of 
the  hyperfine  structure,  it  has  been  an  important  tool  to  determine  the 
magnetic  moments  of  many  nuclei. 

The  analysis  of  molecular  spectra,  however,  has  been  hampered  by  the 
fact  that  the  numerical  computations  necessary  to  obtain,  the  required  accur- 
acy are  beyond  the  limits  of  desk  computing  in  all  but  the  very  simplest  cases 
(such  as  H„) .  These  computations  constitute,  therefore,  a  promising  field 
in  which  to  use  electronic  digital  computers.  They  are,  however,  by  no  means 
straightforward.  The  proper  selection  of  atomic  orbits  requires  a  great  deal 
of  experience  and  thorough  knowledge  in  this  special  fields  careful  and  tedious 
mathematical  analysis  is  necessary  to  keep  the  numerical  work  down  to  what  is 
tolerable  even  for  a  fast  electronic  computer  and  to  guarantee  sufficient  ac- 
curacy by  avoiding  numerically  instable  procedures. 

E.  C  Sahni  and  J.  W.  Cooley  have  attacksd  this  problem,  using  our 
computer  for  the  numerical  work.  While  computer  operation  was  supported  by 
Contract  T)a-36-0'ih-ORI>-l6k6 ,   for  which  this  final  report  has  been  prepared, 
all  analytic  and  coding  work  has  been  supported  by  Contract  NAw  61+75  (between 
the  National  Advisory  Committee  on  Aeronautics  and  New  York  University) 
under  the  terms  of  which  technical  reports  will  be  prepared.  Therefore,  we 
shall  not  report  on  the  results  of  this  work,  but  merely  mention  some  aspects 
of  the  computations  which  may  be  of  more  general  interest. 
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22.51  Decomposition  of  Overlap  Integrals 

All  molecular  integrals  have  integrands  which  depend  on  the  distances 
r^  and  r^  from  f(fo  nuclei  A  and  B.  A  typical  example  is  the  overlap  integral: 

The  meaning  of  the  coodinat«a  can  beat  be  shown  by  a  sketch 


A  ^  (j^  B 


and  in  the  simplest  case,  the  functions  -!>  and  li/,  may  be  a  pair  of  SLATER 
orbitala  such  as  : 


v„  =  i^At  '-" 


e 

'4  ^^  k"*^ 


^.tc 


The  k's  of  "[y    and  Ij/,  are  sometirnes  the  same,  but  more  often  they  are 

different. 

It  would  be  possible,  of  course,  to  carry  out  the  integrations  "by 
brute  force",  i-e.  by  computing  the  integrand  explicitely  in  three-dimension- 
al space  for  a  sufficient  number  of  points  to  achieve  the  desired  accuracy. 
This  method  would  be  neither  attractive  nor  economical.  Fortunately,  the 
angle  CO   is  common  to  the  A  and  the  B  systems  and  the  integration  with  respect 
to  CD  can  be  easily  carried  out.  We  are  thus  left  with  a  two-dimensional 
integral.  For  the  SLATER  orbitals,  this  can  be  further  reduced  to  linear 
combinations  of  one -dimensional  integrals   ^  A^v  ^^   '^n.   ,3^  f<^  exa^rvple  : 


1)  cf.  B.  S.  Mullikan,  C-  A.  Rieke,  D-  Orloff  and  H.  Orloff,  Journal  of  Chem. 

Phys.,  17,  1214-8-1267,  (19i*-9)  where  numerous  further  references  are  given. 
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Sr'Va.s.Ybis)=g^(i^t)^^(i-tf''^[A3Bo-A^lE,-A>S2  +  AoB3] 


k^.i-\<b 


This  reduction  is  achieved  by  introducing  spheroidal  coordinates 


22.52  Computation  of  the  Auxiliary  Integrals  A  and  B 

The  "maater  formulae"  expressing  the  SLATER  atomic  orbit  overlap 
integrals  in  terms  of  A  's  and  B  's  are  so  simple  that  they  can  be  handled 
on  desk  computers.   It  is  worth  while,  therefore,  to  produce  tables  of  the 
A  and  B  integrals.  Preliminary  tables  for  n  =  0  (1)  5  and  a  wide  range  of 
the  argument  have  been  computed  and  more  extensive  tables  are  pl8.nned  for 
publication.  The  following  procedures  have  been  employed  in  order  to  get 
at  least  eight  digit  accuracy: 

(i)  The  A  (p)  can  be  easily  computed  by  partial  integration,  which 
yields  the  recurrence  formula 

This  procedure  is  numerically  stable  since  all  terms  are  positive, 
(ii)  By  partial  integration  of  the  B  integral  one  obtains 
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The  recurrence  formula  contains  both  positive  and  negative  terms  and  is 
unstable  except  for  large  values  of  tke  argument  q. 

(iii)  For  small  values  of  q,  B  (q)  may  be  computed  from  the  power 
series 

PC  0  • 

These  series  are  numerically  stable  but  require  much  more  work  than  the 
recurrence  formulae,  especially  if  q  ia  not  small.  Experinents  have  shown 
that  the  recurrence  formula  yields  the  same  results,  within  the  desired 
eight  digits  accuracy,  if  q>l  and  n<  5,  and, therefore,  it  will  be  used 
for  the  final  tables  within  this  range. 

22.53  Numerical  Integration 

We  have  seen  that  the  overlap  integrals  can  be  expressed  by  the  A^'s 

and  B  's,  which  are  relatively  easy  to  evaluate.  Two  other  types  of  molecular 

n  ' 
integrals,  viz.  the  Coulomb  and  Hybrid  integrals,  can  be  computed  with  the 

help  of  the  auxiliary  function 


where 


^rv(t,T)  = 


1)  Barnett,  M.  P.  and  Coulaon,  C  A-  Phil.  Trans.  Roy.  Soc.  London  2^3  (1951) 
p.  221 
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The  computation  of  the  Bessel  functions  I 


difficulty.  The  recurrence  formula 

^n-l/2 


n+1/2  ^""^  ^1+1/2  Presents  no 


^n+3/2(^^    =  K_, /^(z)    + 


2n+l 


1-1/2 


(2) 


is  convenient,  and  it  ia  stable  for  all  positive  values  of  z.  The  corres- 


ponding recurrence  formula  for  I 

I 

the  power  series 


hl/2 


I   ,  =-(2z) 
n+1/2 


. n+1/2 


(z)  is  not  stable  for  small  z.  Hence 


2L 
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was  uaed  for  z  <  10  and  n  >  1. 


The  integral  Z  p  must  then  be  evaluated  by  numerical  integration.  The 
integrand  ia  positive  throughout,  except  for  t=o,  where  it  vanishes  like 

z     ;  it  will  further  decrease  exponentially  for  x  >oo    .     For  t=  r 

the  integrand  is  continuous,  but  the  first  derivative  is  not.  Thus  the 
integrand  has,  roughly,  the  following  shape 


^  t 


A   number  of  integration  procedures,  such  as  Weddle's  rule,  Gaussion  quadrature 
and  Euler'a  formula,  were  considered,  and  the  last  one  was  finally  chosen,  as 
it  seemed  to  be  not  only  the  simplest  but  also  the  fastest  procedure  for  a  given 
accxiracy.  As  is  well  known,  Euler's  formula  can  be  described  as  integrating, 
for  each  interval,  a  third  degree  polynomial  whose  values  and  first  derivatives 
at  both  ends  agree  with  those  of  the  actual  curve,  while  the  trapezoidal  does 
not  adjust  the  derivative  at  all.  Yet,  the  only  additional  terms  for  Euler's 
formula  are  the  first  derivatives  at  the  end  points  and,  of  course,  at  such 
points  where  the  first  derivative  is  discontinuous.  This  was  easy  to  code  and 
the  results  were  fully  satisfactory,  as  could  be  seen  by  recomputing  the  most 
critical  cases  with  narrower  intervals. 
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22.5^  Coding  Procedure 


As  usual,  we  first  studied  the  question  of  whether  to  use  direct 
machine  language  or  the  FLINT  interpretive  routine.  Since  it  was  evident 
that  serious  scaling  problems  would  arise  and  that  card  punching  would 
take  up  an  appreciable  fraction  of  the  total  machine  time,  it  was  decided 
to  use  FLINT,  but  to  speed  up  the  output  by  changing  the  output  routine. 

The  flexibility  of  FLINT  permitted  a  good  deal  of  experimenting  in 
developing  methods,  and  debugging  was  fast  and  easy.  The  numerical  integra- 
tion, however,  required  quite  a  lot  of  machine  time.  The  tvo  irinermost 
loops  which,  of  course,  take  up  most  of  the  running  tiraft,  vere  then  rewritten 
in  "direct  floating  point",  i.e.  retaining  the  floatiag  pc5.nt  arithmetics 
but  saving  the  time  for  interpretation  and  some  unnessary  no.rmalization  by 
recoding  these  loops  in  machine  language.  This  proved  ■,0  cu^^  machine  time 
down  to  about  one  fourth  without  requiring  too  much  coding  time. 

The  normal  FLINT  output  routine  puts  only  two  floating  point  numbers 
on  each  card  and  punches  only  one  card  at  a  time.  The  new  routine  increased 
these  figures  to  three  and  eight,  respectively.  This  not  only  reduced  the 
output  time  by  appr.  6o^,  but  also  produced  a  format  far  more  suitable  for 
our  tables. 
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22.60  NUMERICAL  EVALUATION  OP  SOME  DOUBLE  AND  TRIPLE  INTEGRALS  ARISING 
FROM  MESON  THEORY 

Originator;  Hiroshl  Suura 

Analysis  and  Coding?  Hans  J.  Maehly  and  Hiroahi  Suura 

Recent  experiments  on  high-energy  electron  scattering  by  hydrogen 
show  that  the  proton  has  an  extended  structure,  as  expected  from  fneson- 
theory,  and  that  the  measured  cross -section  can  be  fitted  approximately 
by  assuming  a  certain  distribution  of  the  proton  c^ge,  and  the  magnetic 
moment  for  this  distribution.  Dr.  Suura  has  derived  i'igoroua  expressions 
for  these  distributions  from  meson  theory.  Each  distribution  is  expressed 
as  the  sum  of  a  double  and  a  triple  integral. 

As  an  example  we  show  one  of  the  triple  integrals 

T/^)  ^_L_  ^  ((fX cr(2;  r_j 1 1 


00 


=  J(M/;U)^H-^Jl-^^-^/^ 


and  the  region  of  integration  is  defined  by  the  inequalities; 

(  q-x  <  y  <  q+x  ^ 

0<x<b;J  i50<2<t" 

I   X  <  y  <  b   J 

It  can  easily  be  seen  from  these  inequalities  that  0  <  q  <  2b. 

The  integrals  were  approximated  by  sums  of  the  integrand  at  points: 

x=  (i.|)cr1  s=l 

y  =     (   j   +  -)'<J    >        i>   J-   Is  being  integer  and  subject  to  inequalities 

?  (        analogous  to  those  for  x,  y,   z. 

z  =     (   k  +|)-c^  J 
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For  convenience,  the  same  Interval  J  was  used  for  q,  i.e.  the  integrals 
were  computed  for 

q  =  (t  +  ^)'(^  ,i  integer ,   0  <  1  <  2n-l 

The  integer  n  and  hence  the  interval  a    vere  not  fixed  constants  but  part 
of  the  input;  the  code  was  run  for  various  values  of  n  and  it  was  found 
that  the  accuracy  obtained  with  cT  =  0.2  was  fully  satisfactory.  The 
computations  were  then  carried  out  for  b  =  5oU  and  b  =  6<,0. 

At  a  first  glance,  the  number  of  points  for  which  ..he  integrand 

k 
must  be  computed  seems  to  be  of  the  order  of  n  ,  requiring  about  one  hour 

computing  time.  However,  the  integral  given  above  as  a  typical  example 

can  be  written  as 

and  CL=  Ji-^k"^'  ^   b=\|TT^  ;  C^nJT^-Z^ 

This  decomposition  brings  the  number  of  integration  poiaxa  down  to  the 
order  of  n  and  the  computing  time  to  3  minutes,  excluding  input  and 
output.  This  is  so  short  that  we  decided  to  use  the  Floating  Point 
Interpretive  Routine  (FLINT)  with  an  estimated  running  time  of  30  minutes. 

The  times  for  coding,  debugging  and  production  are  fairly  tsrpical 
for  FLINT.  Coding,  from  a  finished  flow  chart,  took  one  afternoon.  De- 
bugging was  done  for  n=3,  using  the  FLINT-TRACER.  A  few  trivial  errors 
were  readily  found  and  corrected,  all  this  on  the  evening  of  the  same  day; 
machine  times  approximately  20  minutes.  Production  took  place  the  next 
evening  requiring  another  ^-0  minutes  of  machine  time.  Dr.  Suura  watched 
the  coding  for  the  first  two  integrals  and  then  easily  coded  the  remaining 
two  without  assistance  from  our  staff. 


where 
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23.10  TRAFFIC  SIMJLATIOK  WITH  A  DIGITAL  COMPUTER 
Originator  and  Coder;  S.  Yo  Wong 

Every  one  of  us  has  had  the  experience  that  the  capacity  of  a 
traffic  system  is  determined  by  its  weakest  point.  Millions  may  be 
spent  on  the  construction  of  a  third  lane,  while  the  real  cause  of 
traffic  jams  lies  in  a  quite  restricted  area,  say,  one  traffic  circle. 
It  is  possible,  of  course,  to  find  and  remedy  these  bottlenecks  by  the 
method  of  trial  and  error.  But  it  would  be  much  more  economical  to  re- 
place at  least  part  of  this  experimentation  by  automatic  computing.  The 
two  main  difficulties  hampering  a  computational  approach  are 

(i)  to  determine  the  "input  parameters",  i.e.  the  expected  flow  of 
traffic  at  all  points  and  under  various  conditions  (weekdays,  holidays,  etc.) 

(ii)  to  develop  a  suitable  "model"  of  the  system  under  study  in 
terras  of  discrete  variables  including  the  mechanical  properties  of  the  cars 
and  the  psychological  behavior  and  reactions  of  the  drivers. 

It  has  been  the  purpose  of  the  present  study  to  demonstrate  that  a 
fast  automatic  digital  computer  is  well  suited  to  carry  out  the  actual  computa- 
tion after  the  above  difficulties  have  been  overcome.  A  very  simple  example 
for  which  the  difficulties  (1)  and  (ii)  are  reduced  to  a  minimum  was  carried 
out  numerically. 

23.11  Description  of  the  Simulation  System 

The  read  Section  under  study  was  a  six  lane  one-way  street,  divided 
into  two  three-lane  sections.  The  dividing  island  had  one  gap  which  allow- 
ed for  changing  from  the  right-hand  to  the  left-hand  side.  Each  lane  was, 
for  computational  purposes,  divided  into  unit  blocks  (UB)  of  "unit  car  length" 
(UCL)  3  l8'l4-"  ,'  while  the  tima  step  was  chosen  to  be  one  second.  This  yielded 
a  speed  unit  of  l8'l^"/sec  =  12.5aph.  For  each  time  step  the  care  were  mov©^ 
by  a  discrete  number  of  UCL's,  beginning  with  the  cars  closest  to  the  "exit" 
im  order  to  avoid  double  occupancy  of  UB's.  The  movement  was  further  control- 
led by  some  traffic  regulations  (such  as  v  <  50mph)  and  a 
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certain  probability  distribution  for  the  generation  of  random  numbers 
governing  the  input  of  cars,  their  respective  speed  and  intention  to 
switch  lanes. 

23.12  Conclusion 

The  results  show  that  the  flow  of  traffic  resembles  the  observed 
behavior  to  a  fair  degree,  in  spite  of  all  simplifications,  and  that  the 
distribution  of  transit  times  does  not  depend  too  much  on  the  specific 
set  of  random  numbers  even  for  rather  short  runs  (several  minutes).  Much 
longer  runs  would  be  required,  however,  for  a  realistic,  i.e.  more  complic- 
ated and  extended  road  system.  It  is  obvious,  therefore,  that  faster 
machines  with  larger  high-speed  memories  would  be  needed  to  solve  actual 
problems.  Such  machines  are,  or  will  soon  be  available.  It  will  be  well 
worth  while,  therefore,  to  further  pursue  this  application  of  fast  au- 
tomatic computing  and  to  focus  some  research  on  the  problems  (i)  and  (ii) 
mentioned  in  23-10  above. 
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23-20  A  MIXING  PROBLEM 

Originator  and  Analyst;  H.  H.  Goldstine 
Coder:  Mrs.  M.  Lamb 

This  problem  is  concerned  with  two  compressible  fluids,  one  on  top 
of  the  other  in  a  two  dimensional  container  with  rigid  walls.  The  only- 
exterior  force  acting  is  that  of  gravity.  At  the  beginning  of  the  problem, 
each  fluid  will  be  in  hydrostatic  equilibrium  except  for  a  small  sinusoidal 
perturbation  of  the  interface. 

If  the  upper  fluid  is  lighter  than  the  lower  fluid  at  the  same  pressure 
the  problem  is  almost  trivial.  The  aim  of  our  computation  was  to  determine 
the  sequence  of  events  in  the  reversed  case.  We  shall  expect,  of  course, 
that  the  lighter  fluid  will  "break  through"  at  a  relatively  small  front  and 
then  spread  out  on  the  surface  above  the  upper  fluid. 

The  problem  was  formulated  in  Lagrangian  form.  Accordingly,  at  each 
meshpoint  five  quantities  were  stored,  namely  two  apace  coordinates,  two 
velocities,  and  the  energy.  It  was  found  necessary  to  use  about  600  meshpointa, 
hence  to  store  around  3000  quantities.  As  the  new  drum  was  not  yet  completed 
at  that  time  it  was  necessary  to  store  two  numbers  (of  20  binary  digits)  in 
each  word. 

In  order  to  avoid  the  formation  of  shocks,  which  are  not  relevant  to 
the  main  study,  the  method  of  von  Neumann  and  Bichtoyer  ^was  adapted  to  ©ur 
case. 

It  was  planned  from  the  beginning  to  proceed  step  by  step  from  a  very 
simple  teat  case  (with  the  lighter  fluid  on  top)  to  more  realistic  and  more 
complicated  problems.  However,  great  difficulties  were  encountered  in  form- 
ulating this  problem  and  the  work  was  greatly  hampered  by  the  limited  size 
and  poor  reliability  of  the  old  magnetic  drum.  The  problem  was  finally  abandon- 
ed when  both  Dr.  Goldstine  and  Mrs.  Lamb  left  this  project. 


1)  A  method  for  the  Numerical  Calculation  of  Hydrodynamic  Shocks,  J.  Appl. 
Phys.  21,  232-237  (1950) 
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23 "30  NUMERICAL  UTTEGRATION  OF  THE  NAVIER -STOKES  EQUATIONS  FOR 
COMPRESSIBLE  FLUIDS       ~~ 

Originator,  Analyst  and  Coder:  S.  E.  Lam 

It  is  attempted  to  obtain  a  numerical  solution  to  the  complete 
system  of  partial  differential  equations  governing  the  steady  flow 
field  of  a  compressible,  viscous  and  heat  conducting  gas  in  two  dimens- 
ions. The  physical  problem  selected  for  study  is  the  motion  of  a  thin 
flat  plate  moving  with  hypersonic  speed  at  zero  angle  of  attack.  The 
major  interest  of  this  problem  is  the  behavior  of  the  solution  about  aad 
ahead  of  the  leading  edge,  where  no  analytical  solution,  exact  or  approx- 
imate, is  known.  The  asymptotic  behavior  of  the  solution  at  downstream 
infinity,  however,  is  known  from  the  approximate  boundary  layer  theory. 

The  plate  shall  be  insulated  and  of  zero  thickness,  and  the  "non- 
slip"  condition  on  the  plate  is  assumed  to  be  valid.  The  plate  is  semi- 
infinite. 

23.31  Nature  of  the  Solution 

Since  this  is  a  physical  problem,  the  general  character  of  the 
solution  may  be  inferred  from  the  experimental  data.  We  thus  expect 
that  the  solution  fields  will  be  essentially  undisturbed  far  ahead  and 
far  away  from  the  plate.  A  curved  detached  shock  wave  will  stand  at 
some  distance  in  front  of  the  leading  edge  of  the  plate.  A  boundary 
layer  surrounds  the  immediate  vicinity  of  the  flat  plate,  on  the  surface 
of  which  the  gas  velocity  and  the  normal  temperature  gradient  are  zero. 
At  distances  far  downstream  of  the  leading  edge  the  flow  field  will  be- 
have according  to  the  known  boundary  layer  theory.  Also  the  flow  field 
in  front  of  the  curved  shock  wave  will  be  practically  undisturbed.  The 
information  to  be  obtained  from  the  numerical  solution  concerns  the  shape 
of  the  shock  wave,  the  detachment  distance  of  the  shock  wave  from  the 
leading  edge,  the  flow  field  about  the  leading  edge,  and  the  pressure 
distribution  on  the  plate  surface. 
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23 '32  Equatlona  Governing  the  Flow  Field 

Let  the  plate  be  situated  on  the  positive  X  axis,  with  the  leading 
edge  at  the  origin.  Then 


-I- 


L  ^c*^  >   >ou)  '  -J    Vox-    c-)/    -i  V  <?<    °t)^   J 


vhere        u  velocity  in  +x  direction 

V  "        "  +y 

t  time 

T  temperature 

p  pressure 

^  density, 

y  kinematic  viscosity 

k  conductivity 


and 


C        heat  capacity 
P 

;  =  c  (T)  ^   k  =  k(T)    ;     V  =  V  (T)  ^    ?  =  P  (P'"^^ 


are  known  relations. 
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In  essence  we  have,  therefore,  a  set  of  four  non-lineaj  aimultaneoua 
partial  differential  equations  for  the  four  dependent  variables,  u,  v,  O  ,  T. 
The  character  of  the  equations  is  elliptic  in  space  and  parabolic  in  time,  and 
we  must  impose  appropriate  boundary  conditions  on  a  closed  curve  surrounding 
the  region  of  interest.  The  above  set  of  equations  contains  three  independent 
variables,  x,  y,  and  the  time,  t.  Although  we  are  only  interested  in  the  steady 
behavior  of  the  flow  field,  we  have  been  using  the  unsteady  equations  to  facilit- 
ate a  scheme  of  successive  approximation  to  be  described  later. 

23.33  Main  Difficulties 

As  was  mentioned  above,  we  expect  a  curved  shock  wave  in  the  solutions. 
A  shock  wave  occupies  a  narrow  region  in  which  the  gradients  of  all  dependent 
variables  are  very  large.  If  the  mesh  of  the  flow  field  were  made  fine  enough 
to  permit  an  adequate  representation  of  this  narrow  region,  then  the  number  of 
mesh  points  required  would  be  prohibitively  large.  However,  the  structure  of 
the  shock  is  of  no  interest  in  this  problem.   It  is  only  required  that  the  mag- 
nitudes of  the  dependent  variables  before  and  after  the  shock  can  be  related 
correctly  as  a  function  of  the  orientation  of  the  sho^i;.  Actually,  the  exact 
relations  of  the  dependent  variables  before  and  after  a  shock  as  a  function  of 
the  shock  orientation  are  known  in  analytic  form.  However,  since  the  shape  and 
position  of  the  shock  are  not  known  to  us,  --  indeed,  they  are  the  objective  of 
this  investigation  --we  cannot  take  advantage  of  these  relations. 

Our  first  problem,  now,  is  to  find  a  difference  scheme  which  represents 
the  derivatives  in  the  equations  such  that  the  resulting  system  of  difference 
equations  is  capable  of  yielding  correct  relations  of  the  dependent  variables 
before  and  after  the  shock  for  any  shock  orientation  with  a  very  coarse  mesh, 
perhaps  with  only  two  or  three  mesh  points.  Furthermore,  the  difference  equa- 
tions must  be  stable  in  flow  regions  away  from  the  shock. 

The  second  difficulty  is  the  correct  choice  of  the  initial  and  boundary 
conditions.   If  the  mesh  size  selected  is  fine  enough  to  be  stable,  yet  coarse 
enough  so  that,  with  a  reasonable  number  of  mesh  points,  the  top,  front  and 
downstream  sides  of  the  rectangle  can  be  effectively  considered  to  be  y  =  +oC, 
X  =  -cxD,  and  y:  =  +  00 ,   respectively,  with  the  origin  at  the  leading  edge,  then 
the  boundary  conditions  in  the  steady  state  would  be  essentially  those  described 
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in  23«31«  However,  since  we  choooe  to  approach  th©  steady  state  solution 
from  an  initial  condition  through  a  transient  calculation,  the  appropriate 
boundary  conditions  to  be  imposed  during  transient  are  not  known.  This  pro- 
blem is  still  being  studied  analytically. 

23-3^  Simplified  Problem 

In  order  to  gain  experience  and  insight  concerning  our  system  of  dif- 
ferential equations,  a  much  simpler  problem  was  proposed  and  coded  as  a  trial 
problem.  The  dependent  variable  T,  was  assumed  to  be  constant  throughout 
the  flow  field  at  all  times,  thus  reducing  the  number  of  dependent  variables 
to  three.  Physically  this  corresponds  to  isothermal  flow,  or  a  gas  of  infinite 
heat  conductivity.  This  eliminates  the  possibility  of  a  shock  wave,  therefore 
the  problem  is  much  simpler. 

The  flow  field  is  represented  by  a  32x32  grid,  with  the  flat  plate 
situated  at  the  bottom  of  the  rectangle.  The  leading  edge  is  situated  at  the 
middle,  and  the  flow  is  coming  from  the  left. 

The  computation  was  started  with  the  fictitious  initial  state  of  a  uni- 
form flow  defined  by  u  =  l/U,  v  =  0  and  p  =  \/h   for  the  entire  region. 

Using  the  time -dependent  equations,  one  can  calculate  the  flow  field  at 

t  =  n+1,  knowing  the  flow  field  at  t  =  n.  For  the  initial  condition  of  uniform 

flow  the  gas  velocity  on  the  plate  is  non-zero.  As  the  computation  begins,  the 

boundary  condition  of  non-slip,  i.e.  u  =  0  on  the  plate  surface,  is  gradually 

imposed.  We  used  u   ,  =  u  -l/32,  with  u  =  lA  and  u  =  0  for  n>  8.  The 
n+ln''      o'      n        - 

values  of  u  =  l/U,  p  =  lA,  v  =  0  are  held  rigidly  constant  with  time  on  the 
upstream  and  top  sides  of  the  rectangle.  On  the  portion  of  the  bottom  side 
of  the  rectangle  in  front  of  the  leading  edge,  we  require  v  =  0  and  -^  =  v~  =  0* 
The  purpose  of  the  trial  problem  is  to  study  the  stability  of  our  difference 
equations  and  to  obtain  information  about  the  appropriate  boundary  condition 
for  the  downstream  side  of  the  rectangle. 

If  n  is  the  time  index  and  k  and  m  are  the  x  and  y  space  indices,  res- 
pectively, then  the  difference  schemes  to  represent  the  partial  derivatives  used 
are  " 
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where  V  stands  for  any  dependent  variable.  We  always  chose  idx  =  ^y,  but 
ZLx/z^t  was  left  free  for  adjustment  as  required  for  stability. 

There  is  still  very  little  known  about  the  beundary  condition  t«  be 
applied  on  the  downstream  side,  and  before  we  know  the  analytical  asymptotic 
solution  during  the  transient,  any  inaccuracies  in  the  boundary  condition 
will  generate  errors  which  shall  propagate  upstream.   In  the  trial  runs 
d  V /<3x,  -0  was  tentatively  used  as  a  first  approximation. 

Ay". 

The  result  of  the  trial  runs  for    2  =  2  showed  signs  of  instability 
in  the  difference  scheme  used.  After  about  six  or  eight  time  cycles,  it  was 
found  that  scattered  negative  values  of  u,  v,  and  p  began  to  show  up  in  the  flow 
field,  probably  indicating  numerical  "overflew"  in  the  machine,  i.e.  physically 
unreasonably  large  values.  Up  to  this  time  cycle  the  errors  from  the  approximate 
downstream  boundary  condition  propagate  about  three  mesh  points  upstream. 

23*35  Flans  for  Further  Investigations 

It  is  planned  to  continue  the  search  for  a  stable  difference  scheme  for 
©ur  system  of  partial  differential  equations.  We  hope  to  succeed  in  eliminating 
the  source  of  trouble  at  the  downstream  boundary  by  finding  and  using  the  analjrti- 
cal  as3nnpt©tic  solution  for  x  — >  +  oo.  The  present  knowledge  of  the  appropri- 
ate representation  of  partial  derivatives  in  this  non-linear  equation  system  is 
very  scanty  and  we  have  been  forced  to  use  trial  and  error  methods.   If  a  stable 
difference  scheme  can  be  found,  the  computation  of  solutions  to  this  problem  should 
be  fairly  straightforward.  The  experience  and  insight  gained  in  solving  this 
problem  should  help  us  to  attack  successfully  a  variety  of  compressible  viscous 
flow  problems  in  the  near  future. 
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23.*^0  AUTOMATIC  NETWORK  ANALYSIS 

Originator:  S»  Y.  Wong 
Analyst  and  Coder;  M.  Kochen 

The  analysis  of  complicated  electrical  networks,  for  direct  or 
for  alternating  current,  has  long  been  considered  a  stronghold  of 
analogue  computers  since  the  construction  of  such  network  analysers 
is  pretty  much  straightfoward.  The  application  of  digital  computers, 
however,  offers  the  great  advantage  that  changes  in  the  parameters 
or  even  in  the  connections  can  be  made  much  faster  and  easier,  once  a 
general  code  has  been  prepared.  It  was  the  aim  of  the  authors  to 

investigate  the  practicality  and  limitations  of  digital  computing 

*) 
in  this  field.  The  method  is  explained  in  a  Technical  Beport.  ' 

We  therefore  limit  ourselves  to  a  necessary  correction.. 

As  can  be  seen  from  page  9,  "the  "power  dissipation"  P  is  a  quadra- 
tic function  of  the  voltages  w^  The  Taylor  series  on  page  11  terminates, 
therefore,  with  the  quadratic  term  and  so  far  no  term  has  been  neglected. 

The  statement  on  page  12  "Because  the  Taylor  series  was  truncated 
after  two  terms,  the  above  expression  is  only  a  first  approximation  to 
^V  "  i-B,   therefore,  wrong;  actuallv  ^V   was  deliberately  chosen  to  be 
parallel  to  grad  P. 


*)  M.  Kochen  and  S.  T.  Wong,  Technical  Beport  No.  55-02,  August  1955* 
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23.50     A  TABLE  FCB  CUMULATIVE  BINOMIAL  PROBABILITIES 

Originator:     J.  W.  Tukey 
Coder  :     H.  F.  Trotter 

The  probability  that  in  n  independent  trials   of  an  event  having 
probability  p  of  "success"  the  number  of  successes  will  not  exceed  k  is 
the   cumulative  binomial  probability 

k 
P  (p,  n,  k)   =2Z      (?)   P^Cl-P)"""^ 

This  function  is  not  well  suited  to  ordinary  interpolation  procedures 
and  since  it  has  three  arguments  extensive  tables  are  necessary.  The 
object  of  our  computation  was  to  produce  a  short  table  in  a  modified  form 
such  that  simple  interpolation  will  give  results  of  reasonable  accuracy. 

Define  K(p,  n,  k) ,  the  "equivalent  normal  deviate",  by  the  equation 


J_.(''^  -.<% 


Then 


N(p,  n,  k)  =  2  I  \l   (k+1)  (1-p)  -v/(n-k).p"  } 


is  a  fair  approximation  for  K(p,  n,  k)  and  is  used  as  an  auxiliary  variable 
for  the  table  listing   values  of  K-K  as  a  function  of  p,  E  =  np^ and  N,  with 
the  respective  ranges: 

p  =  0  (0.1)  0.5 

N  =  -1+  (0.5)   -t^ 

2h/-E   =  1   (1)   2k 

As  a  result  of  our  transformation,  this  rather  coarse  mesh  will  permit 
fairly  accurate  interpolation  for  most  values  of  p,  n,  k,  with  the  exception 
of  a  small  region  where  a  finer  mesh  would  be  desirable. 


23.60 

23 '60  Experiments  in  the  Uae  of  FLINT 

During  the  months  of  October,  November,  and  December  in  both  1955 
and  1956  experiments  were  carried  out  to  test  the  effectiveness  of  the 
FLINT  language  as  a  tool  for  the  inexperienced  coder.  Graduate  and  un- 
dergraduate engineers  from  Princeton  University,  without  prior  experience 
in  digital  computing,  were  encouraged  to  code  two  standard  problems  in 
order  to  find  out  what  difficulties  they  encountered.  They  coded  the 
solution  of  n  simultaneous  linear  algebraic  equations  and  the  intergration 
of  an  n   order  ordinary  differential  equation  (Milne,  initial  condition). 
These  routines  are  now  available  as  service  routines. 

The  experiments  showed  remarkably  fewer  coding  errors  than  have  been 
observed  in  machine  language  codes  for  other  machines  —  apparently  attribut- 
able to  a  lack  of  exceptions  in  the  FLINT  language.  The  errors  that  did 
occur  suggest  desirable  changes  for  any  future  external  language  that  may 
be  devised. 
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{ 
23-70  HtSTOEICAL  EPHEMEEIS  FOR  THE  YEAE  -600  to  0 

Originator:  Otto  Neugebauer  and  A-  Sachs 
Analysist  and  Coder;  Bryant  Tuckerman 

Among  the  tools  of  historical  research  ia  the  possibility  of  est- 
ablishing chronology  by  means  of  ancient  astronomical  records.  This 
technique  applies  whenever  it  can  be  shovm,  from  modern  astronomical 
theory,  that  an  observation  or  configuration  recorded  in  an  ancient 
text  could  only  have  oc cured  on  a  particular  date.  This  dates  the 
observation,  and  hence  generally  the  text,  the  calligraphy^ and  any 
contemporary  historical  and  social  records  in  the  text. 

Some  ancient  observations  have  even  been  precise  enough  to  furnish 
improvements  in  the  values  of  some  astronomical  quantities  -  the  empirical 
secular  accelerations  -  which  cannot  be  predicted  from  theory. 

Heretofore  methods  have  not  been  available  for  directly  choosing 
possible  dates  from  the  given  observations.   Instead  it  has  been  necess- 
ary to  choose  provisional  dates,  compute  the  positions  of  the  desired 
bodies  from  existing  theories,  compare  the  results  with  the  recorded  observat- 
ions, and  use  the  results  of  the  comparison  to  estimate  a  new  provisional 
date,  proceeding  in  this  cycle  until  a  satisfactory  fit  is  obtained. 

Each  actual  computation  of  a  planetary  position  has  been  a  fairly 
laborious  process.  A  book  by  P.  V-  Neugebauer,  later  slightly  amended 
by  some  improved  coefficients,  organized  the  procedures  into  a  systematic 
set  of  tables.  While  these  tables  were  a  real  advance  they  still  left  a 
large  amount  of  calculation  to  be  done  in  each  dating  problem,  and  the 
computation  of  a  single  position  might  still  take  about  twenty  minutes 
for  a  person  acquainted  with  the  process. 

There  are  a  large  number  of  astronomical  records  and  tables  surviv- 
ing from  the  region  of  Babylon  in  the  period  from  -600-0.  Professors  0- 
Neugsb&uer  and  A.  Sachs,  faced  with  the  problem  of  dating  these  records, 
proposed  to  Dr.  Goldstine  of  the  Electronic  Computer  Project  the  production 
of  an  Ephemeris  --  that  is,  a  table  of  planetary,  lunar  and  solar  positions 
at  regular  intervals  --  for  the  period  in  question.  With  the  ephemeris 
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available,  the  dating  of  a  record  would  be  reduced  to  lacating  the  desired 
configuration  in  the  ephemeris  by  inspection. 

This  problem  was  thought  appropriate  for  the  Electronic  Computer 
Project  for  several  reasons.  Primarily,  it  was  an  example  of  a  novel 
and  useful  application;  while  the  amount  of  computation  involved  was  of 
such  magnitude  as  to  be  completely  impractical  without  an  electronic  computer, 
it  was  of  considerable  but  not  excessive  size  on  our  computer-  Secondly, 
it  presented  a  number  of  interesting  problems  in  high-precision  calculations, 
error  analysis,  information  handling,  and  methods  of  checking.  Finally,  the 
results  were  of  a  non-commercial  nature  and  would  be  of  permanent  scientific 
value . 

A  thorough  analysis  of  the  problem  was  made  from  original  sources. 
It  was  decided  to  use  five-day  intervals  for  Mercury,  Venus,  and  the  Moon, 
and  ten-day  intervals  for  Mars,  Jupiter,  Saturn,  and  the  Sun.  Except  for 
the  Moon  the  latitudes  and  longitudes  are  given  in  units  of  .01  ,  and  the 
intervals  chosen  are  fine  enough  to  permit  interpolation  with  adequate  accuracy 
by  means  of  Everett's  formula  (for  which  a  simplified  method  of  use  will  be 
furnished).  The  positions  of  the  Moon,  whose  motion  is  too  rapid  to  permit 
satisfactory  interpolation  if  a  reasonable  time  interval  is  chosen,  will  be 
given  in  units  of  .1  . 

Work  on  this  problem  was  commenced  in  the  summer  of  1955*   It  has 
inspired  and  been  paralleled  by  the  development  of  a  number  of  coding  aids, 
such  as  an  assembly  routine,  a  system  of  computer  operation  using  "control- 
cards",  and  an  assortment  of  utility  sub-routines  (see  Part  I).  The  problem 
is  expected  to  enter  production  in  the  summer  of  1957' 

It  is  anticipated  that  the  ephemeris,  when  produced  and  thoroughly 
checked,  will  be  published  by  a  scientific  society  as  a  monograph  of  about 
300  pages  of  tables  with  accompanying  text  explaining  the  method  of  computa- 
tion. To  complement  this  publication,  a  Technical  Eeport  may  be  prepared  for 
Contract  Nonr-1358-(03) ,  by  which  all  mathematical  analysis,  programming  and 
coding  has  been  supported. 


ADDENDUM: 

The  following  two  Technical  Reports  were  issued  under  Contract  No 
DA-36-03l<~ORD-l6li-6 : 

Technical  Report  No  56  -  01,  January  1956: 

"A  METHOD  FOE  FINDING  THE  GENERAL  SOLUTION  TO  AN  ARBITRARY  NON-SINGULAR 
SYSTEM  OF  LINEAR  EQUATIONS  INVOLVING  n-^/2  MULTIPLICATIONS" 

toy 

J.  Paul  Roth 
and 

Technical  Report  No.  56  -  02,  April  I956 

"ALGEBRAIC  TOPOLOGICAL  METHODS  FOR  THE  SYNTHESIS  OF  SWITCHING  CIRCUITS 
IN  n  VARIABLES" 

J.  Paul  Roth 

This  work  is  not  included  in  our  Final  Report  since  full  disclosure  of 
the  results  was  made  in  the  Technical  Reports  and  no  machine  operation 
was  involved . 
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